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Introduction 


Current helicopter rotor blade designs incorporate fiber composite materials as a 
means of controlling weight, deformation, and vibration (i.e., structural tailoring). 
Although fiber composites are orthotropic in material property classification, they can 
exhibit general anisotropic material behavior to applied loads (extension-bend-twist 
coupling) when the fiber orientations do not coincide with the structural coordinate 
system [1]. Exact analytical solutions of this problem based upon three-dimensional 
elasticity are intractable. Three dimensional finite-element modeling may be applied, 
however, this approach is too costly for the discretization necessary for accurate stress 
and deformation determination. 

Instead one-dimensional beam-type models based upon standard or refined 
theories are used. Standard beam theories, which are derived by extending the 
Bernoulli-Euler theory to include extension-bend-twist coupling effects, have been 
developed for thin-wall rectangular cross-sections and general nonhomogeneous cross- 
sections [2,3]. In addition refined beam theories, which include the effects of transverse 
shear deformation, exist for thin-wall single cell cross-sections. All of the existing one- 
dimensional beam-type models yield only gross structural behavior such as force 
resultants, moments, extension, bending rotations, and twisting angles. Stress 
distributions may be calculated according to kinematic (Bernoulli-Euler and Saint-Venant 
torsion) hypotheses, but they do not reveal the inter-laminar shear states. Moreover, up 
to 48 unique constants must be known Tor each different cross-section. These cross- 
section constants provide important information concerning the beam axial, bending, 
torsion, and shear stiffnesses, as well as coupling stiffnesses. But these constants, which 


i.i 


are dependent upon the general warping of the cross-section, can very significantly with 
changes in cross-section geometry and/or material definition. Thus it becomes very 
difficult for an analyst to assess how minor changes in the cross-section definition will 
effect the overall blade behavior. 

Alternatively, Saint-Venant's elasticity solutions for extension, bending, torsion, 
and flexure of a prismatic beam can be used to analyze the cross-section of advanced 
composite helicopter rotor blades. These elasticity solutions accurately describe the 
displacement, stress (including inter-laminar shear), and general cross-section warping 
distribution for a given applied tip-load condition. Thus, these solutions can be used to 
complement existing one-dimensional beam theories by providing a means for: (1.) 
determining the general three-dimensional warping of the cross-section, (2.) calculation 
of the warping-dependent cross-section constants, and (3.) accurate calculation of the 
stress distribution throughout the blade using the calculated shear and moment 
resultants that can be determined from an appropriate derived one dimensional beam 
theory. This beam theory must be derived using assumptions that insure that the 
kinematic and stress fields are fully compatible with the aforementioned Saint-Venant 
elasticity solutions. 


Research Objectives 

During the contracted period, our research was concentrated into three areafi , 

1.) The development of an accurate and a computationally efficient method for 
predicting the cross-section warping functions in an arbitrary cross-section 
composed of isotropic and/or anisotropic materials. This new method involved 
using a "power series^representation for the in-plane and out-of-plane warping 
functions and then solving for the power series coefficients using variational 
principles. Our work developed a theoretical approach and computational 
procedure for cross-sections composed of isotropic materials (Chapter 2) and 
generally anisotropic materials (Chapter 3). In addition, a separate research effort 
was undertaken to develop the v exact"^cross-section warping functions and fully 
coupled three-dimensional displacement and stress field for a tip-loaded 
cantilevered beam having a solid elliptical cross-section and composed of generally 
anisotropic materials (Chapter 4)7 These exact solutions are extremely useful for 
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validating the computationally predicted warping functions and also for correlation 
of the stress and displacements distributions of a one-dimensional beam theory. 

2.) the development of a general higher-order one-dimensional theory for anisotropic 
beams. This theory is used to study the behavior of beams having an arbitrary 
nonhomogeneous cross-section, where the effects of shear deformation and local 
cross-section deformation are included using the aforementioned cross-section 
warping functions. Thus, it is imperative that the beam theory is derived so that the 
kinematic and stress assumptions are fully compatible with those of the Saint- 
Venant elasticity solutions. Numerical results have proven that this new beam 
theory accurately predicts both the displacement and stress distribution (all six 
components). In our work, we first developed a linear dynamic one-dimensional 
theory for isotropic beams, where we studied the static and free vibrational behavior 
so as to assess the importance of the in-plane and out-of-plane warping functions 
(Chapter 7 5). Our results show that the in-plane functions are required for acquiring 
accurate shear stress distributions, whereas only the out-of-plane warping function 
is required for static and free vibration displacement information. Second, we 
developed a general nonlinear one-dimensional theory for spinning anisotropic 
beams (Chapter 6). Our preliminary results show the importance of including both 
the in-plane and out-of-plane warping functions for determining accurate stress 
information and coupled displacement behavior (i.e. static deformed shapes and 
mode shapes). 

^3.) The development of an analytical model for assessing the extension-bend-twist 
coupling behavior of nonhomogeneous anisotropic beams with initial twist (Chapter 
7). A model was formulated, where the displacement solutions are defined with 
pretwist-dependent functions that represent the extension, bending and torsion and 
pretwist-dependent in-plane and out-of-plane cross-section warping functions. 
Numerical results illustrate the strong extension-torsion coupling behavior in thin- 
wall advanced-composite beams as a function of ply angle, initial twist level, and 
initial twist axis location.^ 


In the remaining six chapters of this report, the three different research areas and 
associated sub-research areas are covered independently including separate 
introductions, theoretical developments, numerical results, and references.. This was 
dona because, first, each of the six topics are very independent in their focu_s;and scope 
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and, second, each of the six chapters is a copy of an extended manuscript of an accepted 
and/or published Journal or Conference article. For example, the work of Chapter 2 has 
been published as a conference article [4] and has been extended and accepted to the 
AIAA Journal [5]. The results of Chapter 3, which has been recently completed, has been 
submitted to the International Journal of Solids and Structures [6]. The work of Chapter 4 
has been submitted for publication by The Journal of Composite Materials [7]. The 
results of Chapter 5 have appeared in condensed form in a conference article [8] and will 
appear in extended form in the Journal of Sound and Vibration [9]. The model of Chapter 
6 is being published as a conference article [10] and will be submitted to a Journal when 
additional numerical results are completed. Finally, the formulation of Chapter 7, was 
presented in a conference research article [11] and an AIAA Journal article [12]. In all 9 
articles, the valuable assistance of R.C. Lake and the financial support of the NASA- 
Langley Research Center was acknowledged. 
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Chapter 2: A Power Series Approach for Isotropic Beams with 

Arbitrary Cross-Sections 

Abstract 

The behavior of a tip-loaded cantilever beam with an arbitrary cross-section is studied 
using Saint-Venant’s semi-inverse method along with a power series solution for the out- 
of-plane flexure and torsion warping functions. For complex cross-sections, the 
calculated power series coefficients represent a “best-fit approximation" to the exact 
warping function. The resulting warping functions are used to determine the cross- 
section properties. A new linear relation is developed for locating the shear center, 
where the twist rate is zero about the line of shear centers. Numerical results are 
presented to verify the approach and second provide section data on NACA four-series 
airfoils not currently found in the literature. 

Introduction 

Closed-form solutions for Saint-Venant's flexure and torsion problems (tip-loaded 
cantilever beam) exist for only a few simple cross-section shapes (ellipse, rectangle, 
equilateral triangle) [1-3]. For general cross-section shapes (i.e., cambered airfoils), the 
cross-section dependent flexure and torsion warping functions cannot be determined 
exactly and thus approximate techniques must be used. One proven approach for 
approximately determining the Saint-Venant torsion [4] and flexure [5] functions involves 
the application of the finite element method. In this approach, the general cross-section 
is discretized into triangular and/or quadrilateral subregions (elements) with out-of-plane 
nodal variables that represent the cross-section warping, where the warping distribution 
is determined by applying the principle of minimum potential energy. While the finite 
element approach is well behaved, it has two shortcomings. First a large number of 
elements are required for complex cross-sections, which leads to a large set of linear 
algebraic equations. Second, the resulting array of calculated nodal values provides 
very little physical insight into the warping definition and one typically resorts to graphical 
finite element post-processing techniques to understand the warping distribution. An 
alternative approach, which has been developed by Mindlin [6] for the solution of Saint- 
Venant's torsion problem (generalized plane strain), involves assuming a double power 
series for the warping function. The power series coefficients are determined by solving 
a set of linear algebraic equations, where the number of equations is equal to the 
number of unknown coefficients. Thus, the problem size is independent of the cross- 
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section complexity, and only dependent on the number of terms in the power series. 

The objective of the current investigation is to study the flexure and torsion 
behavior of a tip-loaded cantilever beam with an arbitrary cross-section, where both the 
flexure and torsion warping functions are expressed as a double power series in terms of 
the cross-section coordinates. The coefficients associated with the power series terms 
are determined by solving a set of variationally derived linear algebraic equations, where 
the number of equations is equal to the number of unknown coefficients. For complex 
cross-sections, the calculated coefficients represent a "best-fit approximation to the exact 
warping function which may be an infinite series of transcendental functions. To aid in 
the evaluation of the power series weighted area integrals, the cross-section is 
discretized into a series of triangular subregions, where the integration in each subregion 
is evaluated exactly using Guassian Qaudrature formulas for triangles [7,8]. The triangle 
aspect ratio is not critical as opposed to the finite element method, since the power series 
is a global cross-section function and not a loca] element function. 

Once the flexure and torsion warping functions are known for a given cross- 
section and material definition (Poisson's ratio), then the resulting three-dimensional 
displacement and stress distributions can be used to: (1.) study the overall beam 
behavior, (2.) determine important beam-type section properties including: the torsion 
constant, shear deformation coefficients, shear center location, and shear correction 
factors (for Timoshenko's beam theory [9,10]), and (3.) develop a one-dimensional beam 
theory [11] that includes cross-section flexure and torsion warping effects and is fully 
compatible with the three-dimensional stress and displacements predictions of Saint- 
Venant. 

The determination of the shear center location has been studied by numerous 
researchers [2,3,5,12-15], where the shear center is commonly defined as " the load point 
where the mean value of the local cross-section twist is zero (i.e., local twist rate about 
the section centroid is zero)". Applying this definition to Saint-Venant’s flexure and 
torsion problems leads to a zero twist rate about the centroidal axis, but any other line 
parallel to the centroidal axis, including the line of shear centers, will have a nonzero 
twist rate. This nonzero twist rate for all lines except the centroidal axis occurs because 
the application of a transverse tip load produces a linearly increasing bending stress 
state normal to the cross-section and straining wjthin the cross-section. This straining 
within the cross-section, which causes the particles of the cross-section to translate and 
rotate into an anticlastic surface, also increases linearly from the beam tip. Thus any line 
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that is composed of cross-section particles that are offset from the centroidal axis will 
undergo linearly varying twist (i.e. constant twist rate with zero twist at the beam-tip) as a 
result of an applied transverse tip load. Recently, in [16,17], an analytical approach was 
developed for locating the shear center in thin plate-like cross-sections, where the 
aforementioned definition was modified to be "the load point where the twist rate about 
the line of shear centers is zero". Although the two definitions appear to be very similar, 
this definition insures that the shear center is coincident with the center of twist. 

Moreover, numerical results for thin triangular cross-sections revealed profound . 
differences (40-67% depending upon Poisson ratio) in the shear center locations. In my 
current paper, the shear center location is determined using both the classical definition 
and the more recent definition [16,17], where a linear relationship between the two 
locations is developed that is valid for any cross-section shape and material definition. 

Numerical results are presented to verify the approach and provide new data for NACA 
four-series airfoils. The sensitivity of the section properties (especially the shear center 
location) with airfoil thickness and camber is of interest to aeroelasticians because of the 
profound effect these parameters have on divergence and flutter speed calculations. The 
current work significantly improves upon the fundamental studies of the shear center 
location for solid airfoils [15,18], where these analyses approximately treated cambered 
airfoils as cubic ovals. 

Theoretical Background 


General Beam Behavior 

We begin by considering a cantilever prismatic beam of length L with an arbitrary 
cross-section of area A composed of a homogeneous, isotropic material. A Cartesian 
coordinate system ( x,y,z ) with corresponding displacement components (u,v,w) is defined 
with the origin at the centroid of the root end and the (x,y) axes coincide with the principal 
axes of the root cross-section. See Fig. 1. The beam is subjected to a force with flexure 
components (F x , Py) that act through the centroid of the tip cross-section (z=L) in the x- 
and y-directions, respectively, and an applied torque ( Mz ), where they satisfy the 
following three equations of stress equilibrium; 

P* = J T xz dA , Py ~ | XyzdA , M z = | [xTyz - yt xz )dA . (I.a-c) 


Furthermore, the body forces are considered negligible so that the in-plane stresses (oxx. 
< 7 yy, Txy) are equal to zero and the normal stress is defined as a linear function, following 
Saint- Venant’s assumptions [1-3], 
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where hex and fyy represent the principal moments of inertia about the x- and y-axes, 
respectively. Introducing these assumptions into the stress equil ibri um equations and 
integrating yields the well-known form [1-3] of the displacement components and shear 
stresses 


Elyy 


- + jj* - y 2 j(L - z)j + ^-{ vxy{L - z)j - Qyz + b^z - bey + &i (3.a) 

= -fyjL z 2 - zl + xly2 . *2 1 l - z)l + vxy{L - z)) + Oxz + b$z + b$x + bz (3 .b) 
Elx) <\2 6 2 [ J J Elyy' 

- f )) • -f))"r(x./>-b i x-bsy + b 3 (3.c) 
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(3. d) 
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where, Eis the Young' s modulus of the mater ial, v is the P oisson's rat io, G is the shear 
modulus that satisfies ( = B(2(1 + v )) ), 0 is the beam twist rate about the centroidal axis, 
y/{x,y) is a function that describes warping out of the cross-section plane, and by be are 
integration constants that are specified by defining the fixity of the beam root. In the 
current development, the geometric boundary condi tions are prescri bed by r estraining 
the translational motion of the root centroid (x=>*=z=0) and requiring that the slope and 
twist of the centroidal axis are zero at the beam root, 


u 


v = 


w s 


du 

dz 


dy _ o 

dz ~ dy dx 


(4.a-/) 
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This approach is identical to that of [2,3] and thus £j| =0 (i = 1*6). 

The determination of the warping function (y) is accomplished by applying the 
principle of minimum potential energy 

6n=5U-5W e = 0 (5. a) 


where ( 8U) is the variation of the strain energy given by 


-n 


SU =| | SOzzCzz + fayzYyz + SxxzYxz CM dz , 


(5 .b) 


since (oxx = oyy = r X y = 0), and 6Wq is the variation of the work of external forces that 
results from the applied tractions on the beam ends; 


5W e = | \°zz\ {ZmL) -<Tzz\ (ZM0) 


8ydA 


= xSxydA +^l ySydA. (5 .c) 

1 yyjA Ixx Ja 


Substituting Eqns. (2,3,5 .b, S.c) into (5.a) and carrying out the integration over the beam- 
length 


sn = 0 = GL I ^<5y) 


\ b s A%-r e -^* 2 - y2 i^ VXy] 

Ja \ 


+ #-w I * 1 + *» ■ -fH v *y) - - * 2 > 

ay lay £/ x / 1 &xx \ 2 W 
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( 6 ) 


’P~ EL \ xSYdA--%-EL\ ySydA . 

e, yy Ja e,xx Ja 

An examination of Eq. (6) reveals that (y) can be expressed as a linear combination of 
three cross-section dependent functions that are proportional to the rates of beam 
curvature and twist 
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Wix,y) = Y 3 MO (7) 

Cfyy Clxx 

where (v'l.Y'fc) represent the shear-dependent warping functions and (^3) is the Saint- 
Venant torsion warping function. 

The shear stress distributions of Eqns. (3.d,e) can be expressed in terms of the 
rates of bending curvature and twist by making use of Eq. (7), 


where 



(8 .a,b) 


(8.C./7) 


The twist rate (8) as a function of the applied loads (P x , Py, Mz) can now be 
determined by substituting Eq. (8) into (l.c), integrating and rearranging 


where 
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(/ =1,3) . 


(9.a) 


(9 -b) 
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Torsional Rigidity and Shear Ce nter Location 

The cross-section torsional rigidity [GJ) is commonly defined as the constant (33) in Eq. 
(9 .a). The remaining two constants (a-j,a2) in Eq« (9-3) can be used to locate the classical 
definition of the shear center (x s ,y$); 'the load point that produces a zero mean value 
cross-section twist' (i.e., zero local twist about section centroid. 0 = 0) [2,3,5, 1 2-15]. This 
location can be determined by applying a general flexural force ( Pxs , Pys) through the 
unknown point (x s ,y s ), recognizing the equivalent centroidal forces and moments are 

P x = Pxs , Py « Pys » M z = PysXs-PxsYs, (10.a*C) 


substituting these resultants into Eq. (9.a), rearranging, noting that (0 = 0), and then: 


Xsm *Z- 
s £/«' 


ys =- 


a 1 


El 


'yy 


(10 .d,e) 


Although this definition leads to a zero twist rate about the line of centroids there will be a 
nonzero twist rate about every other line parallel to the centroidal axis as a result of the 
formation of the anticlastic surface. This is observed by calculating the micromolar twist 
rate as 


d(Dz _ 1 

3 z 2 
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d/yz dyxz 
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Applying a force (Pxs, Pys) through the above shear center definition will produce a zero 
micromolar twist rate about the centroid axis (since x=y=0 =0), but a nonzero micromolar 
twist rate (0b) about the calculated line of shear centers (x s ,y s ) that is equal to: 


da>z 

9z 


x-x s 

x-ys 


0c = - 


xs 


El 


'yy 


( v *) 


+ % i,xs ) • 


El, 


'XX' 


( 12 ) 


This nonzero micromolar twist rate is illustrated in Fig. 2.a, where the deformed root (z=0) 
and tip (z=L) cross-sections of the loaded (Pys) cantilever beam are superimposed. 

An alternate shear center location (x s *,ys‘) can also be determined so that the 
application of an applied flexure force (Pxs, Pys) will produce a zero twist rate about the 
calculated line of shear centers and insure that the shear center is coincident with the 
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center of twist. The current model has advantages over the procedure of [16,17] in that it 
uses the Saint-Venant results directly and can be applied to any arbitrary cross-section, 
not just thin cross-sections. This calculation involves finding the load point in the cross- 
section plane that will produce a twist rate about the centroidal axis that is equal to the 
negative of the anticlastically produced micromolar twist rate, thus 



where these twist rates offset each other to produce a zero twist rate about the line of 
shear centers. See Fig. 2 .b. Substituting Eqns. (13) and (lO.a-c) into (9. a), one can 
determine this shear center location as: 

* s " (El xx +vGJ)' Ys = '{Elyy+vGJ) ’ * 14 ' 3 ’* 

or in terms of the previous (classical) definition by making use of Eq. (10 .d,e) and (G/E = 
1/2(1 +v)): 




Ys 


ys 


(1+v) 2/j 


1+7T* 


XX 


(1+v) 2 lyy 


(14.c,d) 


and the percent difference between the two locations can be expressed as: 

*s-*s _ v J Ys-Ys = v J 

X s (1 + v) 2I XX ’ y s (1 + V) 2 lyy 


(14.©./) 


It is interesting to note that (xs*»/s*) will always be less than or equal to (x^/s), where the 
two locations are equal when either the Poisson's ratio is zero, the shear center is 
coincident with the centroid (a 2 =a 3 = 0 ), or the cross-section is very thick so that {Jlkx) or 
(j/ly y) is effectively zero. For most thin solid cross-sections (including NACA airfoils or 
triangles), J/kx - 4, and thus the percent difference when measured from the centroid will 
generally vary from 40% to 67% depending upon the Poisson ratio (assuming 
0.25<v<0.50), which is agreement with the work of [16,17]. 
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Shear Deformation 


U 

s l 

E3 


An examination of the displacement components of Eq. (3 a-c) reveals that the 
calculated centroidal tip displacements ( u{x=y= 0 ,z=L) = P x L 3 /3E/yy, v(x=y= 0 ,z=L) = 
PyL 3 /3E^x) agree with the strength of materials solution, but the additional displacement 
associated with shear deformation does not appear as a result of our original assumption 
(Eq. 4 .a-f) that the slope of the deformed centroidal axis at the beam root is zero 
(64,65=0). This additional displacement can be included by simply rotating the deformed 
beam so that the slope of the deformed cross-section at the centroid (x=y=0) is coincident 
with the x-y plane, and thus the deformed centroidal axis will have a nonzero slope (see 
[ 3 ] for additional details). The rotation angles are equal to the shear strains (ftz.tyz) at 
the centroid of the beam root (x=y=z= 0 ) and thus by combining Eqns. ( 3 , 7 , 9 .a, 10 ) the 
shear angles defined in terms of the applied forces (Px.Py.^t) are: 
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b 4 m Y*Z(o) 


dx 
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dx 
x-y-Oj 


EL 
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dx 
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( 15 .b) 


where the subscript (0) is introduced to symbolize the evaluation of the function at the 
centroid (x=y=0). The final form of the displacement components including shear 
deformation (from Eq. ( 3 .a-c)) is: 


-%*&-** M + ^ - z) < + + 652 


w = 


EL 


' yy 


4 


Lz ' *2^ + ' b * x ‘ b5y 


( 16 . a) 
( 16 . 6 ) 


( 16 .c) 
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where y t 0, 64 , and bs are defined in Eqns. (7), (9 .a), (15.a), and (15 .b), respectively. 


Kinematics of a Compatible One-Dimensional Theory 

A fully compatible one-dimensional beam theory can be developed using the 
following kinematic relations 

u(x,y,z) « U(z) - yS(z) + y/xfx.y), 

v{x,y,z) = V(z) + xS(z) + y/yUy), (17.a-c) 

w(x,y,z) = W(z) - x<Py(z) + y<M*) + yAx.y), 


where (U,V,W) are ^-dependent displacement functions that act along the x,y, and z 
directions, respectively, (&x,<Py,S), are z-dependent rotations about the x,y, and zaxes, 
respectively, and (y/’x.Y'y.V'z) represent cross-section dependent "residual" (or warping) 
displacements of the beam. The in-plane functions (yx,Yy) which are associated with 
formation of the anticlastic surface, can be neglected by assuming (oxx= 0 yy=*xy=O) and 
thus using a one-dimensional constitutive model: 


Ozz - Eezz - E 


BW d*y d0xy 

dz d z y dz 


T yz = G Kyz ~ G 


da dv 




dz dz 


dy 


(18.a-c) 


fxz= Gy xz = G\-y^ + ^-<Py + ^ 


If one assumes that (U,V,W) and (<P*, <fy t O) represent the displacements and 
rotations about the centroid [9], then the corresponding definition of the one-dimensional 
out-of-plane warping function (from Eq. (16.c)) is 


y z = y{x,fi - b 4 x - bsy . 


(19) 
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u 




where {y/) is given in Eq. (7). 

Alternatively, based upon the work of [10], one can assume that {U, V,W) and 
are the mean displacements and the mean rotations of the cross-section 

U, V, W M i-l (u, v, u) dA , <Kx - -H y wdA , <Py - • -f- xwdA , (20 .a-e) 
a )a Wa 

and the correct form of the out-of-plane warping is obtained by substituting Eq. (16.c) into 
(20.a-e) and then into (17.c): 

y/z m yr ~ -j— I yvdA --*- 1 xydA --]r| yrdA , (20 .f) 

IxxJa lyyjA a Ja 


where (y/j is given in Eq. (7). 

The development of the equations of motion and the corresponding boundary 
conditions using the kinematic relations of Eqns. (17.a-c) can be determined using 
Hamilton's principle, where the bending and torsion related section constants are 
dependent upon the one-dimensional out-of-plane warping function {y z )- The complete 
details of this refined model can be found in the second part of this paper [1 1]. 


Finally, a set of linear equations can be developed that relates the kinematic 
description of shear strain and twist rate to the shear (Px.Py) and torsion (A/fz) resultants, 
which can be used to (1) transform the warping function definition (Eq. (7)) to a 
kinematically scaled function, and (2) provide valuable one-dimensional cross-section 
constants. Substituting Eqns. (18.b,c) into Eqns (l.a-c) and carrying out the integration 
over the cross-section results in 


*11 

*12 

0 

*21 

P 22 

0 

*31 

^32 

1 


Px 

Py 

M z 


S^^ 0 S 13 

0 S22 S23 

0 0 S33 



\ 


dV 
9 z 


+ &x ) . 


90 

Bz 


( 21 ) 
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where the coefficients (fljj) and (Sjj) are defined in the Appendix. Multiplying Eq. (21) by 
the inverse of [f?] results in the following set of linear equations 


’ 1 

Px 


GAk u 

GA ftj2 

-GAy k-|3 

Py 

! = 

GA k ^2 

GA k 22 

GAx k 23 

M z 

1 1 


-GAy Ai 3 

GAx ^23 

GJ 


3 U . 

2 K+* 

dz 


30 

3z 


\ 


/ 


( 22 ) 


where (/rn, A12. ^22) are the shear correction coefficients needed for Timoshenko's beam 
theory [9] and (/C13, /C23) are the shear correction coefficients for coupling between flexure 
and torsion. This approach for flexure-torsion behavior Is an extension of the method 
developed in [10] for uncoupled bending only. 


Solution Procedure 

Warping Function Determination via Power Series 

The warping function (y) of Eq. ( 7 ), which is dependent upon both the cross- 
section shape as well as the material properties, can be determined by solving a set of 
variationally derived algebraic equations based upon the principle of minimum potential 
energy. In this development the warping function is defined as a power series: 


00 00 

Y(x,y) = X Z CmnX m y n - Cqq , ( 23 ) 

m-0 0 

where c^n are the unknown coefficients and the rigid body translation coefficient (cqo) is 
not included since it was accounted for in Eq. (3.c) by (63). If one assumes a finite series, 
then the above equation can be written in matrix form as 

Y(x,/I - [M*.y)](c) (24) 

where [N(x,y)] is an array of the power terms, {c} is an array of unknown coefficients, and 
the array sizes are dependent upon the selected polynomial order. For example, if a 
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cubic polynomial was selected, then based upon Pascal's triangle there are nine terms, 
and the above arrays have the form: 


[M*.y)j = { *. y, x 2 , xy , y 2 , x3, x 2 y, xy 2 , y 3 ) , 
|c | * | C 10 . C oi . c 20 . c 11 , C 02 • c 30 • c 21 • c 12 • c 03 } 


(25. a, b) 


A set of linear algebraic equations for determining the coefficients {c} can be obtained by 
substituting Eq. (24), into Eq. (6) and taking the variation with respect to the unknown 
coefficients (6 y= [/V(x,y)]{8c}): 


K |c) -IFhHF^O) 

where the stiffness matrix is defined as 

d 


M - gl 




(26. a) 


£{N(x,rt] + -^[N(x,yi\ ^A/(x,y)] dA (26. b) 


the force matrices are presented as 
[Fvi » E 


[F c ] = GL 


|jM W. 


f y[N] T dA. 

0 , 

JA 

- 

-vxy •] 

4 

( y2 -* 2 ) 

-sM 

-vxy 


(26. c) 


cM , (26. c/) 


and 



(26. e) 


The coefficients {cO)} associated with the unit warping function (v^i) in Eq. (7) are 
determined by setting {Q} 7 '- {1,0,0}. Similarly, the coefficients {c( 2 )} for (y 2 ) and {c( 3 )} 
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for (y/3) are determined by performing analyses with {Q}^ = {0,1,0} and {Q} T = {0,0,1}, 
respectively. Thus, the complete warping function distribution for the three cases can be 
written matrix form as; 

( Vi . V 2 . V 3 ) - W*’?> ]{(e (1) ) ■ jc (2) } ,{c (3 ») } ■ (27) 


Computer Program 

A computer program was written where, first, the b ou ndary of a general cross- 
section is defined using (n) coordinate points with (nj straight fine segments connecting 
the points. Second, the cross-section is discretized into (n) triangular subregions, where 
one edge of a triangle is a boundary line segment and the other two edges connect the 
end-points of a boundary line segment with the user-defined cross-section origin. Thus 
all of the subregions have one corner that is defined at the origin. See Fig. 3 .b for an 
example of a rectangle defined using four triangular subregions. Third, the cross-section 
centroid and principal axes are calculated and then the cross-section coordinates and 
applied forces are transformed to the cross-section principal axes. Fourth, the area 
integrals (Eqns. 26 .b-d) for each triangle subregion are evaluated using exact Gaussian 
Quadrature formulas [8], where the cross-section power series polynomial can be user- 
defined. Fifth, the complete cross-section stiffness and force matrices are formed by 
simply adding together (not finite element type assembling) all of the triangular subregion 
matrices. Sixth, the coefficients for each of the three cases of {Q} are determined. 
Seventh, the calculated coefficients along with the power series polynomial definition are 
used to determine the shear stress distribution and the cross-section properties (shear 
center location, torsion constant, shear correction factors, etc.). Finally the calculated 
values are transformed from the cross-section principal axes back to the user defined 
coordinate system. 

This approach strongly differs from previous finite element based approaches [4,5] 
in that the global matrix size is defined by the assumed polynomial order and not the 
complexity of the cross-section. Moreover, cross-section cavities can be easily treated by 
simply subtracting off the triangular subregions that define the cavity. The aspect ratio of 
a triangle subregion is not critical, since the power series is a global cross-section 
function and not a local subregion function (i.e., finite element method). 


Numerical Results 

Prismatic cantilever beams with three different cross-section types are studied to 
first verify the current approach (ellipse, rectangle), Fig. 3, and second illustrate important 
results not found in the literature (NACA 4-digit airfoils), Fig. 4. The beam material 
properties are defined as (E=1 .0, v^O.333). 

Ellipse 


Initially the current approach was verified by studying the behavior of cantilever 
beams with elliptical cross-sections having a wide range of aspect ratios (O.OI^JVa^lOO), 
where (a) and ( 6 ) are the semi-axes length in the xand / directions, respectively. See 
Figure 3. a. Each elliptical cross-section was discretized using 90 points on the cross- 
section (i.e. 90 triangular subregions) and the warping functions were defined using a 
cubic polynomial (9 unknown coefficients) with a corresponding exact Gaussian 
Quadrature formula [ 8 ]. In Fig. 5, the variation of the nonzero nondimensionalized 
coefficients for an applied bending curvature rate (P x /Efyy) and twist rate ( 6 ) are 
presented as a function of aspect ratio (b/a), where the circles are the current calculated 
power series solution and the bold solid lines are the exact solution from [2]. The flexure 
solution in the x direction is composed of only three nonzero terms (yr ^ = c-io* + c^xy 2 + 
c>j qx 3 ), where the remaining 6 calculated coefficients are equal to zero. The Cio term is 
proportional to the shear strain at the beam root (tfcz(o)) and thus is used as a measure of 
shear deformation in the x-direction (£ 4 ), from Eq. (15.a): 


cio _ 1 <Vi 

a 2 a 2 dx 

x=y~0 


v*z(o) El yy 

a*P x 


/?4E/yy 

a*P x 


(28. a) 


Applying a force in the /-direction (Py), would produce only three nonzero coefficients; 
coi. C 21 , and C 03 , where again the remaining 6 coefficients are equal to zero. The 
coefficient coi is proportional to (fy Z (o)) and represents a measure of shear deformation 
in the y-direction (£> 5 ): 


flU _ 1 dV '2 
b 2 *" b 2 

x-y* 0 


ryz {o) Elxx ^ bsEl xx 
b 2 Py b 2 Py 


(28 .b) 
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The twist rate dependent warping function is composed of only the bi-linear term (yr 3 = 
ciixy), where the remaining 8 coefficients are equal to zero and cn goes to zero as the 
cross-section becomes a circle (b/a=1). The torsion constant (GJ^a^) was calculated for 
each aspect ratio and was found to be in exact agreement (identical to 8 decimal places) 
with the closed-form solution of [2]: 


GJ * G ic — (29) 

a 2 + to 2 

which was expected since the current torsion warping function is in agreement with the 
published solutions. The shear center was calculated and found to be coincident with 
the centroid for all aspect ratios (x s *=x s =0, ys*=ys=0). Finally, the shear correction factors 
(ki 1 .^ 22 ) were calculated and compared to the results of Cowper [10] (see Table 1). The 
current predictions are in near exact agreement over a broad range of aspect ratios, 
where it is interesting to note that for a force acting through a very thin ellipse (for 
example Py with b/a~ 0) the shear correction factor approaches zero while a force acting 
through a thick ellipse (P x with b/a= 0) the shear correction factor approaches 0.917184. 

Rectangle 


The behavior of cantilever beams having rectangular cross-sections was also 
studied in order to further validate the current approach. A wide range of aspect ratios 
(0.01 <b/a^1 00) were investigated, where each cross-section was discretized using the 
four corner points (four triangular subregions). See Fig. 3 .b. From [2], the exact solution 
of the x-dependent bending curvature rate warping function (y^) is defined as: 


V 1 



sinh(/?xf) 
to- 
ri 3 cosh(nxf) 
D 



cos(rwr£). (30) 

b 


In order to compare the current power series predictions with the above infinite series of 
transcendental functions, a Taylor series expansion was performed on Eq. (30) and the 
first three nonzero terms were found to be: 


Y'l 


a 2 



(1 +v) + 4 vS(-2) 


x- ^ 1 + 4S (0 ) |xy 2 + Jj2 + v + 4vS (0 ) |x 2 + . . . (31 .a) 


TO) 


■il 


'( 0 ) 
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where 



(31 -b) 


The torsion rate warping function (y/3) is also described by an infinite transcendental 
series [2], but as (b/a) approaches either zero or infinity the series reduces to simply 
(V3=cnxy) with cn=-1 or cn=1, respectively. Moreover, the torsion warping function for 
a square cross-section (b/a=1) can be expressed using a Taylor series expansion as: 



The calculated torsion rate warping function (v^) is used to determine the torsion rigidity 

GJ = ag = Gk t (2af(2b) (33) 


where (/ct) is the torsion constant. 

Two cross-section aspect ratios (b/a - 1, 100) were initially studied to assess the 
convergence of the calculated coefficients and the warping-dependent cross-section 
constants as a function of the polynomial order for the warping function. In Tables 2 and 
3, the first three nonzero flexure coefficients, the first torsion coefficient, the shear 
correction factor, and the torsion constant are presented as a function of power series 
order and solution matrix size. In addition, reference values for the calculated values are 
presented, where the three flexure coefficients are determined using Eq. (31), the torsion 
coefficient is taken from Eq. (32), the shear correction factor (fcn) is taken from [10], and 
the torsion constant (/ct) is taken from [2,3]. From these tables, it is obvious that the 
integrated cross-section constants converge to the reference values much quicker than 
the actual power series coefficients. This occurs because the calculated coefficients 
represent a "best-fit" of the user-defined polynomial to the transcendental series, and 
, changing the order of the polynomial will change the magnitude of the calculated 
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coefficients, but it will have virtually no effect on the integrals of these functions. Thus, if 
one is only interested in warping related cross-section constants, then a low-order power 
series polynomial can be used, but if one is interested in the details of the warping 
function, then a much higher-order polynomial is required. 

In Fig. 6, the first three nonzero nondimensionalized power series flexure 
coefficients (symbols) are presented along with the Taylor series representation of Eq. 
(30) given by the bold solid lines. The power series solution produces near exact 
agreement over a broad range of aspect ratios (0.0l£b/a<100), where a ninth-order 
polynomial (54 unknowns) was used for the warping function. As (d/a) approaches 100 
(transversely loaded plate-type cross-section), the power series predictions deviate from 
the Taylor series representation. This can be traced to the fact that the infinite 
transcendental series in Eq. (30) converges slowly for large ( b/a ) and the selected ninth- 
order polynomial can not accurately represent this behavior. Thus, one would need to go 
to an even higher-order polynomial for this severe aspect ratio. The shear center (x s ,y s ) 
and shear correction factors (fcn, /C 22 ) were calculated for the entire range of aspect 
ratios and it was found that the shear center was always located at the centroid and the 
shear correction factors were always equal to 0.85105, which is in exact agreement with 
[10] using v= 0.333. The calculated torsion constant (/q) is presented in Table 4 as a 
function of aspect ratio, where the current results are in near perfect agreement with [2]. 

NACA 4-Dig it Airfoils 

The final set of beam cross-sections that were investigated included six NACA 4-digit 
airfoils of different thickness (NACA 0006, 0012, 0018) and camber (NACA 2512, 4512, 
6512). The numbering system for these airfoils is based upon section geometry [19], 
where the first digit indicates the maximum value of the mean-line ordinate in percent of 
chord (c), the second digit indicates the distance from the leading edge to the maximum 
camber location in tenths of chord, and the last two digits indicatejhe maximum thickness 
(fmax) in percent of chord. Two of the studied airfoils, NACA-0012 and NACA-4512, are 
presented in Fig. 4, where a second coordinate system ( x,y) is introduced with the origin 
taken as the leading edge. Each airfoil was discretized using 95 points on the cross- 
section boundary (i.e. 95 triangular subregions), the warping functions were modeled 
using a 9th order power series polynomial (54 unknown coefficients), and the numerical 
integration was performed exactly using a 52-point Gaussian Qaudrature formula [8]. 

The calculated section properties for the six airfoils are presented in Table 5. The 
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first five parameters represent the chord normalized geometric section constants and the 
sixth parameter (/3) is the rotation angle from (x,y) to the principal axes (x,y) with counter 
clockwise defined as positive. The remaining nine parameters represent the torsion- and 
flexure-dependent values. The torsion coefficient (^)* which is nearly independent of 
airfoil thickness and camber, is found using the calculated torsion constant ( GJ) 

GJ s= 3$ = G k{ (fmaxjpc • (34) 

The chord normalized flexure dependent coefficients (cio. 001 ) ar® a ' s ° presented and 
are used to provide a measure of chord-wise and thickness-wise shear deformation (see 
Eqns. 28.a,b). From the coefficient (cio), it is readily apparent that the application of a 
force in the chordwise direction (P x ) will produce nearly constant shear-deformation 
regardless of airfoil thickness or camber. Whereas the shear deformation associated 
with a thickness-wise force (Py) is highly dependent upon the airfoil thickness but only 
slightly dependent upon camber. The shear center locations are also presented for the 
two definitions using the second coordinate system (x,p). Both definitions locate the 
shear center ahead of the centroid, where the difference between the two definitions is 
nearly 50% (measured from the centroid with v=0.333) in the x-direction and a minimal 
amount in the /-direction (cambered airfoils). It is interesting to note that the shear center 
moves rearward (closer to the centroid) by either increasing the thickness or camber. 

The difference that these two shear center locations have on the torsional divergence 
speed of a straight uniform aircraft wing can be studied using (from [20]): 

U D = jL J a J -r , (35) 

2L Y cea^p/2) 

where (p) is the air density, (ao) is the lift curve slope, and ( e ) is the distance from the 
elastic axis (line of shear centers) to the quarter-chord. The ratio of the divergence speed 
( Uo *) using the corrected shear center location (aligned with the center of twist) based 
upon [16,17] to the divergence speed (Ob) using the classical definition [12], while 
holding everything else constant, is equal to: 



where the results for the six different airfoil sections are given in Table 5. Thus it is 
apparent that the new (corrected) shear center location will predict divergence speeds 
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that are significantly lower (6-12%) than those using the classic shear center definition 
and this reduction is larger for thin airfoils with little or no camber. 

Lastly, the shear correction factors (/fi 1 ,^ 22 ) are presented. Increasing the airfoil 
thickness will produce a minor decrease in (fcn) but increase (£ 22 ) significantly. This 
observation Is in close agreement with the thin elliptical cross-section results (Table 1, 
0.01<fe/a<0.20). Introducing camber will significantly reduce both shear correction 
factors (ki 1 ,^ 22 ). 


Conclusions 

The flexure-torsion behavior of a tip-loaded cantilever beam with an arbitrary cross- 
section is studied using Saint-Venant's semi-inverse method along with a power series 
solution for the out-of-plane flexure and torsion warping functions. The power series 
coefficients are determined by solving a set of variationally derived linear algebraic 
equations. For complex cross-sections, the calculated coefficients represent a "best-fit 
approximation" to the exact warping function. A new linear relation is developed for 
locating the shear center using the Saint-Venant flexure and torsion solutions, where the 
twist rate is zero about the line of shear centers (not the centroidal axis). In addition the 
kinematic relations for a fully compatible one-dimensional beam theory are presented, 
where the calculated current flexure and torsion warping functions are fully integrated 
into the development (see part II [11]). Numerical results are presented for three different 
cross-sections (ellipse, rectangle, NACA four-digit airfoils). For elliptical cross-sections, it 
was shown that the calculated coefficients, as well as all of the section properties, were in 
exact agreement with existing elasticity solutions. For the rectangular cross-section, it 
was shown that the calculated power series coefficients represent a "best-fit" to the 
transcendental functions and a low-order polynomial can be used if only warping-related 
section properties are desired, whereas a higher-order polynomial is required if the 
warping function is to be studied in detail. Finally for NACA four-digit airfoils, the shear 
deformation and shear correction factor associated with a thickness-wise force (Py) Is 
highly dependent upon the airfoil thickness but only slightly dependent upon camber. 

The x-direction shear center location is ahead of the centroid with the difference in the 
two definitions being nearly 50% (for v^0.333) when measured from the airfoil centroid, 
and increasing either the airfoil thickness or the camber will move the shear center closer 
to the centroid. These differences correspond to a 6-12% decrease in the divergence 
1 speed for the corrected shear center definition versus the classic definition. 
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Table 1 : Shear correction factors for elliptical cross-sections (v-0.333). 


b'a 

kll 

kll 111] 

k22 

0.01 

0.917181 

0.917181 

0.004777 

0.10 

0.916854 

0.916854 

0.309770 

0.20 

0.915869 

0.915869 

0.602246 

0.50 

0.909272 

0.909272 

0.829613 

1.00 

0.888864 

0.888864 

0.888864 

2.00 

0.829613 

0.829613 

0.909272 

5.00 

0.602246 

0.602246 

0.915869 

10.0 

0.309770 

0.309970 

0.916854 

20.0 

0.105567 

0.105567 

0.917102 

100.0 

0.004777 

0.004777 

0.917181 
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Table 2: Calculated nonzero flexure and torsion power series coefficients, shear correction factor, and 
torsion constant for a square cross-section (b/a-1.00) with (v*0.333). 


Polynomial 

Order 

Matrix 

Size 

Flexure - x direction 

Torsion 


C30 

c 12 

*11 



2 

5 

■ 



1.000 

■I 


4 

14 


•0.3888 


0.85105 

B3I 


6 

27 

1.2324 


0.1018 

0.85105 

mm 

0.141 

8 

44 

1.2340 

-0.3690 

0.1108 

0.85105 

1 .574 

0.141 

Reference 

1.2340 

m 

0.1114 

0.85105 

1.574 

0.141 


2. Z5 

















Table 3: Calculated nonzero flexure and torsion power series coefficients, shear correction factor, and 
torsion constant for a thin rectangular cross-section (b/a«1 00) with (v^O.333). 


Polynomial 

Order 

Matrix 

Size 

Flexure - x direction 

Torsion 

cio/a 2 

C30 

c 12 

*n 

C11 

fa 

2 

5 

-554.0 

- 


1.000 

1.00 

0.333 

4 

14 

0.7783 

-0.3888 

-0.1663 

0.85105 

1.00 

0.333 

6 

27 

1 .2483 

-0.1118 

-0.1670 

0.85105 

1.00 

0.333 

8 

44 

1.1130 

-0.2334 

-0.1514 

0.85105 

1.00 

0.333 

Reference 

0.9943 

■m 


0.85105 

1.00 

0.333 



















Table 4: Torsion coefficients (kO for rectangular cross-sections (v*0.333). 
bfa kt kt [2] 


1.00 

0.14058 

0.141 

1.20 

0.16613 

0.166 

1.50 

0.19578 

0.196 

2.00 

0.22871 

0.229 

2.50 

0.24940 

0.249 

3.00 

0.26336 

0.263 

4.00 

0.28086 

0.281 

5.00 

0.29137 

0.291 

10.00 

0.31297 

0.312 

oo 

0.33333 

0.333 
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Table 5: Section properties of NACA four-digit airfoils (v*0.333). 



NACA- 

0006 

NACA- 

0012 

NACA- 

0018 

NACA- 

2512 

NACA- 

4512 

NACA- 

6512 

centroid 

x/c: 

0.42067 

0.42067 

0.42067 

0.42061 

0.42039 

0.42006 

y/c: 

0.00000 

0.00000 

0.00000 

0.01520 

0.03037 

0.04550 

A/c 2 

0.04106 

0.08213 

0.12319 

0.08219 

0.08238 

0.08270 

lxx/c 4 (10-5) 
lyy/C 4 ( 1 0' 3 ) 

(3 (degrees) 

0.84944 

2.26500 

0.00000 

6.79550 

4.52990 

0.00000 

22.9350 

6.79490 

0.00000 

6.97050 

4.53690 

0.41351 

7.49940 

4.55840 

0.82850 

8.39180 

4.59370 

1.24740 

kt 

0.15642 

0.15386 

0.14986 

0.15396 

0.15427 

0.15479 

cio/c 2 (10* 1 ) 

1.89102 

1.89302 

1.89633 

1 .89609 

1.90553 

1.92117 

coi/c 2 (10’ 4 ) 

8.15829 

32.5947 

73.1685 

32.5100 

32.2590 

31.8141 

Xs / c : 

0.33935 

0.34093 

0.34350 

0.34399 

0.35233 

0.36402 

rJ 

Ys/c: 

0.00000 

0.00000 

0.00000 

0.02183 

0.04378 

0.06595 

x s /c: 

0.36634 

0.36711 

0.36838 

0.36871 

0.37319 

0.37985 

rJ 

ys/c: 

0.00000 

0.00000 

0.00000 

0.02196 

0.04398 

0.06614 

kn 

0.91545 

0.91462 

0.91323 

0.90181 

0.86340 

0.80324 

k22 

0.07740 

0.24212 

0.40041 

0.22576 

0.19355 

0.16671 

U D */Ud 

0.87640 

0.88120 

0.89350 

0.89100 

0.91140 

0.93710 
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3.) Elliptical and rectangular cross-sections with aspect ratio (b/a). The rectangular 
cross-section reveals the four triangular subregions. 
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6.) Comparison of the first-three nonzero calculated flexure power series coefficients 
with the Taylor series expansion of the results from [2] ( — ) for a rectangular 
cross-section. 
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Chapter 3: A Power Series Approach for Generally Anisotropic 

Beams Arbitrary Cross-Sections 

Abstract 

The behavior of a tip-loaded anisotropic cantilever beam with an arbitrary cross-section 
is studied using Saint-Venant's semi-inverse method along with a power series solution 
for the local in-plane and out-of-plane deformation warping functions. The power series 
coefficients are determined by solving a set of variationally derived linear algebraic 
equations. Using the resulting three-dimensional displacement solutions, the shear 
deformation associated with applied tip loads is investigated as well as the shear center 
location. Both of the extended definitions reveal the linear dependency of the shear 
center location with beam-length. Numerical results are presented for three different 
cross-sections (ellipse, triangle, NACA-0012 airfoil ) and two different materials (Al 6061 - 
T6, off-angle high-strength graphite/epoxy fibers). 

Introduction 

Closed-form solutions for Saint-Venant's problems (tip-loaded cantilever beam) 
exist for only a few simple isotropic homogeneous cross-section shapes (ellipse, 
rectangle, equilateral triangle) (Sokolnikoff, 1956) and one anisotropic homogeneous 
cross-section (ellipse) (le and Kosmatka, 1992). For general cross-section shapes, the 
local deformation functions of the cross-section cannot be determined exactly and thus 
approximate techniques must be used. One proven approach for approximately 
determining these local deformations in isotropic cross-sections (Herrmann, 1965; 

Mason and Herrmann, 1968) and anisotropic cross-sections (Kosmatka and Dong 
(1991)) involves the application of the finite element method. In this approach, the 
general anisotropic cross-section is discretized into triangular and/or quadrilateral 
subregions (elements) with in-plane and out-of-plane nodal variables that represent the 
local in-plane deformations and out-of-plane warping. But the finite element method 
requires a large number of elements for complex cross-sections, which will lead to a 
large set of linear algebraic equations (thousands of unknowns). Moreover, the 
calculated array of nodal deformations provides little insight into the deformation and 
warping distribution over the cross-section and thus one must resort to graphical finite 
element post-processing techniques to understand this distribution. An alternative 
approach, which has been used by Mindlin (1975) for the solution of Saint-Venant's 
isotropic torsion problem and by Kosmatka (1992) for the isotropic flexure problem, 
involves assuming a double power series for the each of the local in-plane deformations 
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and the out-of-plane warping. The power series coefficients are determined by solving a 
set of linear algebraic equations, where the number of equations is equal to the number 
of unknown coefficients. Thus, the problem size is independent of the cross-section 
complexity, and only dependent on the number of terms in the power series. 

The objective of this paper is to develop a method for studying the behavior of tip- 
loaded anisotropic beams with general cross-sections using Saint-Ven ant's semi-inverse 
method, where the local deformations of the cross-section are expressed as a double 
power series in terms of the cross-section coordinates. The coefficients associated with 
the power series terms are determined by solving a set of variationally derived linear 
algebraic equations, where the number of equations is equal to the number of unknown 
coefficients. For complex cross-sections, the calculated coefficients represent a "best-fit 
approximation" to the exact warping function which may be an infinite series of 
transcendental functions. To aid in the evaluation of the power series weighted area 
integrals, the cross-section is discretized into a series of triangular subregions, where the 
integration in each subregion is evaluated exactly using Guassian Qaudrature formulas 
for triangles (Dunavant, 1985). The triangle aspect ratio is not critical, as opposed to the 
finite element method, since the power series is a global cross-section function and not a 
local element function. 

Numerical results are presented for three different cross-sections (ellipse, triangle, 
NACA-0012) and two different materials (Al 6061 -T6, off-angle high-strength 
graphite/epoxy fibers) to first validate the approach, second prove convergence of 
warping related cross-section parameters (torsion constant, shear center location, shear 
deformation), third present important behavioral data not currently found in the literature 
and fourth investigate the sensitivity of the shear center location with cross-section 
shape, beam length, and material definition. 

Theoretical Background 

We begin by considering a cantilevered prismatic beam of length L with an 
arbitrary cross-section of area A composed of a homogeneous, rectilinearty anisotropic 
material. A rectangular Cartesian coordinate system ( x,y,z ) with corresponding 
displacements {u,v,w) is established with the origin at the centroid of the root end and the 
(x,y) axes coincide with the cross-section principal axes. See Fig. 1. The constitutive 
relations for the material are given by: 
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M - [c](e) . 

(«} = [S](<t) . 

[C] = [S ]-' , 


(1 .a-c) 


where [C] and [S] are fully populated symmetric matrices with 21 distinct elements and 
the stress and strain arrays are given as: • 


(crj (<Txx . Gyy . Gzz - * yz . T xz t Txy} . 

(e) = { £ xx . £ yy < £ zz . Yyz . Yxz . Yxy) • 




At the root end, the beam is fully fixed. Within the framework of the Saint-Venant 
problems, this condition cannot be described on a point-wise basis and the equivalent 
statement at the centroid (x=y=z= 0) can be used: 


du _ dv _ dv du 
dz dz dx dy 


At the free end, tractions are applied which reduce to an equivalent force P and moment 
Mwith respect to the cross-section centroid. The force Pand moment M can be 
decomposed into flexure components: Pxand Py, an extensional component; Pz, bend- 
ing moments: Mx & nd My, and a torsion moment; Mz . As a result of the applied tip loads, 
five of the stresses are independent of z and the sixth stress (ozz) has flexure 
components which vary linearly with z 


<*zz 


*yy ixx 


z + ozA x <y) 


(3) 


where /xxand lyy are the area moments of inertia about the x and y axes, respectively, 
and o°zz is associated with extension, bending, and torsion. Introducing these 
assumptions into the stress equilibrium equations and integrating yields the following 
displacement and strain components (see Kosmatka, 1986; Kosmatka and Dong, 1991; 
for further details): 
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u = 


Px z2 


2Elyy\ 3 


{z - 3L) - ^-yz(z - 2 L) + j Vi x2 - Vz y 2 j(z - i)J + 2£jrJ M r + ’ e Y z 


2EI 


XX 


[2v y x + v 6 y)y(z- L) + ^-yz{z- 2L)J + V'xfry) 


(4.a) 


Pj 

2 El, 


y|(z-3L)-fx^-2L) + {v 2 y2 . v^z- + f^ 2 + 0*z 

J( v 6 x + 2 v 2 y)>(z - L) + ^-X2(z - 2L)J + ly/x.y) 


(4.6) 


w = - 


2EL 


^-j j v 5 x + v 4 y )x(z - L) - X2(z - 2L)j - ^-| My + ^-M^Jxz + ^— | M x + ^-M 2 jyz 


■ tifcih* + v 4 y}y(^ - <-) - - 2/-)) - ^{f p* 


+ ^.p y .p z ) Z+ y^( Xi y) 


(4.C) 


and 




e yy 


. .Pg2 iz . L) .^ A2 . L) ^M 

Elyy E! xx oy 


(4.d) 


(4.e) 


Ezz = " 


2E/i 


f-{v 5 x + v 4 y+ 2(L - Z )\x---f-{v 5 x + v 4 y+ 2(L - z))y 

/yy v Ztzlxx 


1 


ryz = • 


yxz = - 


Elxx 
Px |v e 


M x + ^-M z jy -^{% + yM z )x-^(v 5 P x + v 4 P y - 2P Z ) 


JjJ^ x+ V2 y+ v 4 (z- l )) x + 2^{ Vix2 ' V2y2 " 2 v 4>< z ' L )j + 6x + 


I y Vc 

— ^ v 1 x + -^-y + 

c/ xx t 4: 


v 5( z -*-)}> 


2EI: 

Px 


(tyz 

dy 


2EI 


yy 


Vi x 2 - V 2 y 2 + 2vsx(z - L) 




(4 .f) 
(4-g) 
(4-/7; 
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where, A is the cross-section area, 6 is the twist rate of the beam about the centroidal axis 
(z), E (= I/S 33 ) is commonly called Young's modulus, and v| are cross-coupling 


coefficients defined as: 



(4.J) 


Here vi and v2, are the usual Poisson coefficients, and v4, V5, and ve, express the 
three-dimensional extension-shear coupling that can occur in a completely anisotropic 
body. The remaining functions (y/*, y/y, y r£) represent local cross-section (x,y) dependent 
deformations and are unique for each cross-section shape and material configuration, 
and are linearly proportional to the six applied loads and twist rate (see Kosmatka and 
Dong, 1991) 


{Wx . Yy >Yz) - X (vOr(i) » Yy(i) . (4-k) 

/-I v ' 

where Q are the components of 

{Q) T = [Px, Py , Pz, M x , M y , M z , 6) . (4./) 

In the current development, the "unit" local cross-section deformations (yr X (i)» Yy(\)< 
Yz(i)) are assumed to have the form of a power series: 

00 00 

Yx(l) = X X a mn {l) xmyn , 
m = 0 n = 0 

00 00 

Yy{i) ~ X X bmn(f)X m y n > (5 .a-c) 

m=0 n=0 


V'z(z) » X 

m= 0 


X C mn (/) 

n = 0 


x™y n 


3.5 


where (a mn (i),bmn(i).Cmn(i)) are unknown coefficients that depend upon the cross-section 
shape, material properties, and load-type ( Q ), and the subscripts (m) and (n) 
correspond to the order of (x) and (y), respectively. The four rigid body motions of the 
cross-section (three translations, rotation about the z-axisj for each of the seven cases 
are constrained by setting (aoo = boo - $00 = 0) and requiring that (aoi = bio). Assuming 
that the series is finite, Eqns (A.k) and (5.a-c) are combined to form: 


where 




{V) T = [Yx, Yy < Yz }. 


[H] = 


[N(x,y)] 0 0 

0 [N{x,yj] 0 
0 0 [N{x,rf] 


(6.a) 


(6 .b,c) 


and [*y] is comprised of seven columns of unknown coefficients that have the form: 



For example, if a cubic polynomial was selected, then based upon Pascal's triangle, 
[A/(x,y)] has 10 terms: 


(6.d) 


[N{x,y)] = { 1, x, y, x 2 , xy, y 2 , x 3 , x 2 y, xy 2 y 3 J 
and {a (j)}, {b (j)} and {c (j)} for the / th column of [‘f'] has the form: 


(7. a) 



%) (/) * a 10 (/) > ^1 (/) > a 20 (i) • a 11 (/) • ^2 (i) > a 30 (/) - a 21 (/) > a 12 (/) • ^3 (/) j • 
boo (/) . bio (/) . boi (/) , b 2 o (/) . bn (,) , bo 2 (/) . ^30 (/) . *fei w . bi 2 w , bo3 (/) ) .(7-b-d) 
( Coo (/) • c 10 (/) . $01 (;) . $20 (/) . C 11 (/) • $02 (/) . $J0 (/) > $21 (/) - c 12 (/) - $03 (/) } - 
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and the aforementioned three rigid body translations and one rigid body rotation are 
constrained using standard finite procedures after the cross-section model is fully 
assembled. 

The strain array {e} of Eq. (I.e) can be obtained in terms of the matrix of unknown 
coefficients [*P], the applied forces and moments, and the centroidal twist rate by 
substituting Eq. (6.a) into Eqns. (4.c/-/): 


where 



W- 


B 




Q . 


ax 


o 


tea 

ay 


0 

0 

0 


0 


0 


0 


3y 3x 


<Wxyj\ 

ay 

ax 


o 


(M 


(8 .b) 


and [ F c ] is defined in the Appendix. 

The magnitude of the unknown coefficients in [*F] can be determined by applying 
the principle of minimum potential energy: 


sn = 5U - 5W e = 0 


(9.a) 


where 5U is the variation of the strain energy; 


5U 


{ f {&) r [C](c}<M dz, 


(9.6) 


and 8Wq is the variation of the work of external forces that results from the applied 
tractions on the beam ends; 

8W e = J [t xz 8Yx + TyzSyy +C zz 8Yz } | ( ZxLj dA ' J [ t xzSYx + TyzSYy +<Jzz$Yz I ( Z=0) dA 
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which reduces to (Kosmatka and Dong, 1991): 


(9 c) 


5W e 


x 8y z dA +^ 1 | ySyzdA. 
l W)A lxx JA 


(M 


A set of linear algebraic equations for determining the seven "unit" unknown 
deformation coefficients is obtained by substituting (8 .a), (9 .b), and (9 .d) into (9. a), 
integrating over the beam volume, and taking the variation with respect to the unknown 
coefficients: 


where 



(10. a) 


(lO.b) 


0 


0 


0 0 0 0 0 


[Fiv] ■ L 


0 0 0 0 0 


(10 .C) 


— I x[N{x,/)]dA -M y[N{x,y)]dA 0 0 0 0 0 

'yyjA xx Ja 


and 


[Fc] 


= l f[e] T tc|l 


F c dA . 


(10. d) 


with [Fc] also being defined in the Appendix. The final form of the local cross-section 
deformations (Eq. (4./c) ) are determined by solving (Eq. (lO.a) ) for pF] and substituting 
the results into Eq. (6.a). Similarly, the stress components (Eq. (1 .d)) can be written in 
terms of (Q), using Eqns. (1 .a) and (8.a), as; 


M = [o][Q]' 

as 


(12. a) 


where 



(12 .b) 


The centroidal twist rate 6 (O7) can now be determined in terms of the three 
applied forces and moments by substituting the fourth (r yz ) and fifth (r xz ) rows of Eq. 
(12.a) into the cross-section torsion moment equilibrium equation 

xXyz-yr xz dA , (13. a) 



integrating, and rearranging to get: 

6=3 - | P x + ^ 2 Py + 3%Pz + 3^M X + 3$My + 3qM z , 


where 


and 


3k = 



a 6 


1 

37 ' 


3k 



x(c 4 k) - y(o5/c)j dA 


(k = 1-5) 


(/c=1,7) . 


(13 .b) 


(13.c) 
(13 .d) 

(13.e) 


The (/) and (j) subscripts of <r,j in Eq. (13.e) correspond to the row and column position in 
[o] (Eq. (12.b)). The coefficients ai-a6 are all independent of the beam-length because 
from the original assumptions, the shear stress distribution is only a function of the cross- 
section coordinates (x) and (y). Moreover, the torsion stiffness is commonly defined as: 
GJ= 1/a6 (see Kosmatka and Dong, 1991). 

Lastly, the local deformations, the strain array, and the stress array can be 
expressed in terms of only the three applied forces and moments by combining Eq. (13.b) 
with Eqns. (6.a), (8.a), (1 1) and (12.a): 


* -mm- 


M - 6 r + r c r ° • 


(14.6) 


(ct) = C \B *P + F C T 


(14.c) 


m- 


(14-cO 


a 2 a 3 a 4 a 5 % 


of = {P* Py, Pz, M x , My, M z ) . 


(14.e) 


Behavior of an Anisotropic Beam 

The general behavior of a cantilever anisotropic beam having an arbitrary cross 
section can now be studied using the displacement (4.a-c) and stress distributions (14.b) 
along with the calculated twist rate (13.b) and the cross-section deformations (14.a). In a 
previous paper, (Kosmatka and Dong, 1991), a detailed discussion was presented 
covering the extension, bending, torsion, and flexure behavior of anisotropic cantilever 
prismatic beams based upon Saint-Venant solutions. In the current paper, we will focus 
our discussions on two topics: shear deformation and further issues concerning the shear 
center location. 


An examination of the transverse displacements (u,v), from Eqns. (4 .a,b), reveals 
that applying either a bending moment or a flexure force (P x ,Py) will produce 

centroidal tip components (x=y= 0, z=L) that agree with the standard (isotropic) strength of 
material solutions 
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(15.a, t») 




u{x=y=0,z=L) 


Px L 3 My L 2 
3 Elyy 2 Elyy’ 


v{x=y=0,z=L) 


PyL? My L 2 
3 E/yy ' 2 E/yy’ 


W 


6 





P 


PI 

pis? 


and applying an extension force (P 2 ) will produce only an axial component (w). But the 
transverse components associated with shear deformation are not included in Eqns. (4 .a- 
c) because the fixed root boundary was defined by setting the deformed centroidal axis 
slope to zero {du/dz-dv/dz= 0). These additional transverse components can be included 
by simply rotating the deformed beam so that the slope of the deformed cross-section at 
the centroid (x=y=0) is coincident with the x-y plane, and thus the deformed centroidal 
axis will have a nonzero slope at the origin. These rotation angles are equal to the shear 
strains (ftz.tyz) at the centroid of the beam root and can be found evaluating the 4 th and 
5 th equations of Eq. (14.b) at (x=y=z=0). For example the rotation angle, about the y-axis, 
associated with shear deformation in the x-z plane is equal (from Eqns. (4 .h), (5.c), and 
(14.a))to: 


YxZ(o) - 


Yxz 


x=y=Q 


dx 
x=y = 0 


6 

X 

/= 1 


C 10(/) +a i c 10(7) 


Qi 


(16) 


Similarly the rotation angle, about the -x-axis, associated with shear deformation in the y- 
z plane is equal to: 


dy z 

Yy z io) = — 


3y 
x=y = 0 


2 |^)1 (/) +a /P01(7)P/ 
/-I 


(17) 


The final form of the displacement components including shear deformation are 

P -*-lzi{z - 31.) - |Ly 2 (z - 2L) + j v, *2 - j(z - t)J + - JpjlWy t }z2 - 6yz 


u = - 


2 Elyy\ 3 


2 El 


xx 


{2v-)X + v 6 y}y(z - L) + ^-yz(z - 2L)j + y XZ(0) z + y/x(x,y) 
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(18. a) 


. v ’ •i^lf {2 ' 3i - ) 'f ,(2(z ' 2t)+ { V2y2 ' Vix2 l (r ' i) )' zki Mx *f M ¥* exz 

- 2^j( VeX + 2 v 2 y)i z • L ) + - 2L )j + fy z (o) z + vofoy) 

(18.6) 

^ - -^{^^Ay}^-L)-^z-2L)).^M y + fM^ Z + ^M x + fM 2 )yz 

■ 2^!i vsx + v>y ^ z ■ L) ■ yi ~ z ' 2i) ) ■ Mf Px + f Py ' p *) z 

- Yxz {0) * - ryz {0) y + Yz(x,y) 

(18 .C) 

where 0, yx. Yy> Yz . ttz( 0 ). and tyz( 0 ) are defined in Eqns. (13.6), (14.a), (16), and (17), 
respectively. 

Shear Center Location. Line of Shear Centers 


For prismatic cantilever beams that exhibit less than generally anisotropic 
behavior (V 4 ,V 5 = 0 ), the shear center is a property of the cross-section and independent 
of beam-length (line of shear centers is parallel to the centroidal axis). For this class of 
materials, a classic definition has been presented for locating the shear center (Griffith 
and Taylor, 1917) as ' the load point that produces a zero mean value cross-section twist 
rate (i.e., zero twist rate about the centroidal axis/. Attempting to extend this definition to 
a beam composed of a generally anisotropic material (V 4 ,V 5 * 0 ) leads to a shear center 
location that is a function of the cross-section shape, material definition, and is linearly 
dependent upon beam-length (Kosmatka and Dong, 1991). Thus, the line of shear 
centers for a generally anisotropic beam is straight, but it is not parallel to the centroidal 
axis. 


We can study this phenomena by calculating the micromolar twist rate for a 
particle in the beam: 


dcoz _ 1_ 
dz 2 


d/yz dyxz 


dx dy 


= e - 


2 EL 


yy 


j V 6 x + 2v;>y + v 4 (z - L)} + 25^{ 2viX + v 6 y + v d z ■ L )) 


( 19 ) 
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w 


m 




w 


where the ( z ) dependent terms are associated with anisotropic 'bend-twist' coupling as a 
result of applied flexure. The micromolar twist rate about the centroidal axis (x=y=0), to 
be consistent with the work of (Griffith and Taylor, 1917), is: 

— 'Mz-L) + £gkz-L). (20) 

Next, we apply a general tip flexural force ( Pxs , Pys ) through the unknown shear center 
location (xs.ys). where the equivalent centroidal forces and moments are defined as 

Px - Pxs . Py = Pys . Mz = Pys x s~ PxsYs ■ (21 -3-c) 

Substituting Eqns. (21. a-c) and (13.b) into Eq. (20) results in 

77 = 1 • <Ws - 2^ U - ■ «) + Py\a 2 + Ws + - «) • (22) 

Since the micromolar twist rate (dcoz/dz) varies linearly with (z) for a generally anisotropic 
beam, it is not possible to locate the shear center so that the twist rate about the 
centroidal axis is zero independent of beam axial position. Thus the Griffith and Taylor 
definition can not be implemented if (V4,V5*0). Instead, we recognize that since the twist 
‘ rate varies linearly with (z), then the micromolar twist (o>z) will vary quadratically with (z). 
Thus, the best that can be achieved is to "locate the shear center (x s ,y s ) so that there is 
zero micromolar twist at the beam root (z=0) and zero twist in the cross-section plane that 
contains the applied load (for a tip load, z=L)". This is accomplished by requiring that 



mm 


Substituting Eq. (22) into Eq. (23), carrying out the integration, and solving, produces the 
shear center location for the beam tip cross-section ( z=L ) that is independent of the 
magnitude of the applied loads: 


x s (z=L) = x sL = 


a 2 


. Ysjz 
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, y s {z=L) = y sL = ± (a 1 + 
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The shear center location at the root of a generally anisotropic beam (or a very short 
anisotropic beam (L= 0)) is equal to the classic Griffith/Taylor definition for an isotropic 
cross-section (x so = -az!a§, y so = a^/ae). The shear center location in the beam-tip cross- 
section (Eq. (24.a,b)) is composed of two terms; one which is independent of the beam- 
length (i.e. the classic Griffith/Taylor definition) and one which is linearly proportional to 
the beam-length, the material properties (V 4 ,V 5 ), and the ratio of the torsion stiffness to 
the bending stiffness {V{aeEl xx ), 1/(a 6 E/ yy )). This development and discussion is 
consistent with the shear center location proposed by le and Kosmatka (1992) for a 
generally anisotropic beam with an elliptical cross-section. 


An alternate shear center definition has been proposed by Reissner (1989, 1991) 
using thin-plate theory so that "the application of an applied flexure force will produce 
zero twist about the calculated line of shear centers f and thus insure that the shear center 
is coincident with the center of twist. Recently, Kosmatka (1992) applied this definition to 
Saint-Venant's flexure and torsion problems and developed a linear relationship for 
calculating this new shear center location in an isotropic beam with an arbitrary cross- 
section: 


*s = 


_§ 2 _ 


ag + 


El. 


xx 


y s = 


a i 


ag + 


E/, 


yy 


{25. a, b) 




where the difference between this definition and the classic Griffith/Taylor definition for an 
isotropic cross-section is an additional term in the denominator. 


Now to apply this definition to a generally anisotropic beam with an arbitrary cross- 
section, we again apply a general flexure force ( Pxs , Pys) through an unknown shear 
center location (x s l*. y S L*) in the tip (z=L) cross-section, where the equivalent centroidal 
forces and moments are given in Eq. (21. a-c) and the micromolar twist rate about the 
unknown shear center location is determined by substituting Eqns. (21. a-c) and (13.b) 
into Eq. (19): 


do) z = p 


dz 


xa 


a i - *6 y S L ■ 2ElJ\ V6X ' sL + 2v2y * L + V ^ Z ' L )) 


(26) 


+ P. 


ft 


a 2 + a 6 x sL + 2 £^j 2 v i *sL + WsL + v 5{z * /-)} 
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Once again, the micromolar twist rate (dcoz/dz) varies linearly with (z), and thus the best 
that can be attained for an anisotropic beam is to locate the shear center so that the 
micromolar twist (coz) is zero at the beam root (z=0) and beam tip (load plane, z=L) and 
varies quadratically between the root and tip. Substituting Eq. (26) into Eq. (23) and 
carrying out the integration over the beam-length, results in two coupled linear algebraic 
equations that are used to solve for the new shear center location in the beam-tip cross- 
section: 
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(27) 


An examination of the above equation reveals that: (1.) this new shear center location for 
anisotropic beams will identically reduce to the results of Eq. (25) for isotropic materials 
(V4,V5,V6=0), (2) the form of the length-dependency effects is identical to that of Eq. (24), 
.which was developed by extending the Griffith/Taylor definition for anisotropy, and (3.) 
the presence of (v 6 ) introduces coupling between the (x) and (y) locations. The coupling 
associated with (V6) is unique in that it is not present in the extension of the Griffith/Taylor 
definition (Eq. (24)) and furthermore, (v6) type material coupling can not be included by 
studying plate-type theories, which make use of plane stress assumptions. Finally, the 
shear center location (xso*, y$o*) in the beam root cross-section is determined by solving 
Eq. (27), where ( L ) is set to zero. An example of the two different 'line of shear centers' 
definitions are presented in Fig. 2, where it is possible that the shear center in the tip 
cross-section plane can be well outside of the cross-section planform. 


Computer Program 

A computer program was written where, first, the boundary of a general cross- 
section is defined using ( n ) coordinate points with (n) straight line segments connecting 
the points. Second, the cross-section is discretized into ( n ) triangular subregions, where 
one edge of a triangle is a boundary line segment and the other two edges connect the 
end-points of a boundary line segment with the user-defined cross-section origin. Thus 
all of the subregions have one corner that is defined at the origin. Third, the cross- 
section centroid and principal axes are calculated and then the user-defined cross- 
section coordinates are transformed to the cross-section principal axes. Next, the area 
integrals (Eqns. 10 .b-c/) for each triangle subregion are evaluated using exact Gaussian 


Quadrature formulas (Dunavant, 1985), where the cross-section power series polynomial 
can be user-defined. Fifth, the complete cross-section stiffness and force matrices are 
formed by simply adding together (not finite element type assembling) all of the triangular 
subregion matrices and the three rigid body translations and rigid body rotation are 
constrained. Sixth, the coefficients for each of the seven cases of {Q} are determined. 
Seventh, the calculated coefficients along with the power series polynomial definition are 
used to determine the shear stress distribution, the constants (a { (/= 1-7)), the cross- 
section properties (shear center location, torsion constant, shear deformation, etc.), and 
transform the seven sets of calculated power series coefficients to six sets associated 
with the six applied loads. Finally, the calculated values are transformed from the 
principal axes back to the user defined coordinate system. 

This approach strongly differs from our previous finite element based approach 
(Kosmatka and Dong, 1991) in that the global matrix size is defined by the assumed 
polynomial order and not the complexity of the cross-section. Moreover, cross-section 
cavities can be easily treated by simply subtracting off the triangular subregions that 
define the cavity. The aspect ratio of a triangle subregion is not critical, since the power 
series is a global cross-section function and not a local subregion function (i.e., finite 
element method). 

Numerical Results 


Prismatic cantilever beams with three different cross-section types are studied to 
validate the current approach (ellipse), prove convergence (triangles), and illustrate 
interesting beam behavior not found in the literature (triangle, NACA-0012 airfoil). See 
Fig. 3.a-c. Two different materials are considered including; an isotropic material (Al 
6061 -T6, E=69 GPa, v^O.333) and a transversely isotropic material (unidirectional high- 
strength graphite/epoxy fibers, (Table 1)) with G 23 = ^22/(2(1+V23))- Generally 
orthotropic or anisotropic beam behavior is introduced by orienting the material reference 
frame (1,2,3) associated with the graphite/epoxy fibers relative to the beam coordinate 
frame (x,y,z) using (a) and (/3), where the angles are defined in Fig. 4 and the resulting 
transformation between the two orthonormal coordinate systems is given as: 



sin(a) cos(/J) 
cos(a) cos(/j) 
- sin(0) 


sin(a) sin(/3) 
cos(a) sin(/J) 
cos(/3) 


cos(a) 
- sin(a) 
0 



( 28 ) 
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The resulting 21 distinct flexibility coefficients (Sjj) are determined using standard 
techniques (see Lekhnitskii, 1963) and are presented in detail in a recent paper (le and 
Kosmatka, 1992; Appendix B). Aligning the fiber axes with the beam coordinate axes 
with (a=/J=0) will result in transversely isotropic beam behavior with V4=V5=V6=0. 
Rotating the material axes about the (y) axis (set /3=0, and vary a), will produce 
orthotropic beam behavior with (v5*0) and v4=ve=0. Similarly, rotating the material 
axes about the (x) axis (set /3=90°, and vary a), will produce (v4*0, v5=vs=0). Finally, 
rotating the material fibers in the cross-section (x-y) plane (set a=90°, and vary P) results 
in (v4=v5=0, V6*0). 


Elliptical Cross-Section 

The current approach was initially verified by calculating the local cross-section 
deformations associated with applied flexure (P y ) for anisotropic cantilever beams with 
elliptical cross-sections having three different aspect ratios ( bla = 0.1, 1.0, 10.0), see Fig. 
3.a, and comparing with the exact results presented by le and Kosmatka (1992). 
Anisotropic behavior was introduced by rotating the graphite/epoxy fibers in the x-z plane 
(0°<a<90°, P=0°) so that vi, V2, and vs are nonzero. Each elliptical cross-section was 
discretized using 90 points along the perimeter (i.e. divide the cross-section into 90 
triangular subregions) and the local deformations of the cross-section (y/x, Yy, Yz) were 
modeled using cubic polynomials (Eq. 7 .a-d). Thus the resulting matrix equation of (Eq. 
lO.a) had 26 unknowns. The calculated power series coefficients were found to be 
comprised of only nine nonzero values that have the form 

Yx = P>|aoiy+ a 2 ix 2 y+ ao 3 y 3 ) . 

y f y = Pyjbio* + £>12 *y 2 + £>30* 3 J , (29. a- c) 

Yz = Pyjqn/ + C2^*y+ C03Y 3 ) . 
where these nine coefficients are expressed from Eq. (14) as: 
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a m n - a m n ( 2 ) + a 2 a m n ( 7 ) 


bmn — {bm n{2) + a 2 bm n(7)) » (29. d-f) 

Crnn = | c /77 n ( 2 ) + a 2 c mn ( 7 ) j - 

and the subscripts (2) and (7) are associated with the second and seventh column of [¥']. 
The remaining 17 coefficients were always identically equal to zero independent of 
aspect ratio and orientation angle. In Figs. 5-7, the variation of the nine nonzero local 
deformation coefficients are presented as a function of orientation angle (a) for the three 
aspect ratios, where the the exact solutions (le and Kosmatka, (1992)) are represented 
by solid lines and the current approach is represented with circles. From these figures it 
is clear that the current approach can reproduce the exact local deformation results over 
a broad range of aspect ratios and orientation angles. Using these local cross-section 
deformation functions, the torsion constant, shear deformation, and shear center location 
were calculated and also found to be in complete agreement with the exact solutions. 
Readers interested in further discussionsjnvolving the variation of cross-section warping, 
stress distribution, and section properties with fiber orientation angle and elliptical cross- 
section aspect ratio are referred to (le and Kosmatka (1992)). 

Triangular Cross-Section 

A second set of homogeneous isotropic and anisotropic cantilever beams having 
triangular cross-sections were analyzed to first prove convergence of the cross-section 
parameters with increasing power series polynomial order and second illustrate 
interesting section property information not found in the literature. The triangle 
represents an interesting cross-section shape because even though it is geometrically 
simple (3 corner points, three triangular subregions), closed-form torsion and flexure 
solutions for the local cross-section deformations exist for only the isotropic equilateral 
triangle, whereas the local cross-section deformations for any other aspect ratio ( b/a ) are 
represented by an infinite series of transcendental functions. For these cross-sections, 
the current approach represents a "best-fit" to the infinite series, where almost all of the 
calculated coefficents will be nonzero. As the order of the power series polynomial is 
increased, the calculated coefficients may vary slightly, but the calculated cross-section 
integrals (section properties) will experience virtually no change. In this study, a 9 th 
order power-series polynomial was used for each of the three local deformation functions 
(161 total unknown coefficients) and the numerical integration was performed exactly 
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using a 52-point Gaussian Quadrature formula (Dunavant, 1985). A second coordinate 
system (x,y) is introduced in the triangle (Fig. 3 .b), where the origin is located at the mid- 
length of the base (b) and (x) bisects the triangle. 

To study convergence of the section properties with polynomial order, we initially 
analyze an isotropic (Al 6061 -T6) cantilever beam with a thin triangular cross-section 
{Ua= 10, b/a= 0.1). In Table 2, the key section parameters are presented as a function of 
polynomial order and matrix size ( [K] ). The normalized shear center locations, x s /a and 
%*/a, are presented in the ( x,y ) system, where the two values are significantly different for 
this thin cross-section, but both approaches exhibit monotonic convergence. The torsion 
constant (GJ = 1 / 33 ) is also presented, where again the solutions converge quickly. 

Lastly, the ratios of the centroidal tip displacement associated with shear deformation to 
the total centroidal tip displacement for applied flexure loads (P x ) and (P y ) are presented, 
where the magnitudes of the ratios are given, from Eq. (18.a-c), as: 

^shear _ ^ z (o)^ ( v shear _ . (30. a, b) 

Utotal , p /_ 3 ''total p 1 3 

* z <°> L + 3 Ei^y *W- + 3 ££ 

As expected, shear deformation is a much larger effect for flexure loads applied in 
x-direction because the effective beam-length aspect ratio is much shorter in the x-z 
plane (L/a =10) compared to the y-z plane (L/b = 100). Both values exhibit monotonic 
convergence, where it is interesting to note that the x-direction value converges quickly 
using a low polynomial order (3). A second convergence study was performed using the 
same geometric beam dimensions, but generally anisotropic behavior was introduced by 
orientating unidirectional high-strength graphite/epoxy fibers with a=/3= 30°. The 
calculated section properties are presented as a function of polynomial order in Table 3. 
The beam root and tip shear center locations are presented using both approaches, 
where it is observed that: (1.) the location converges monotonically with increasing 
polynomial order, (2.) the root shear center locations are within the cross-section, (3.) the 
tip shear center locations are well outside of the cross-section, and (4.) and the result 
obtained by extending Reissner's approach is closer to the centroid (xs/a=0.333), more 
conservative, than the result obtained by extending the Griffith/Taylor approach. The 
torsion constant (GJ) and the ratios of the centroidal tip displacement associated with 
shear deformation to the total centroidal tip displacement for applied flexure loads (P x ) 
and (Py) are also presented, where again these parameters converge monotonically. 


3.19 


In addition to the convergence study, three studies are presented that investigate 
the variation of the shear center location with cross-section aspect ratio and material 
properties. In Fig. 8, the sensitivity of the shear center location with aspect ratio (b/a) and 
Poisson's ratio (v) is presented for an isotropic cantilever beam. The bold solid line and 
the thin solid lines represent the shear center locations based upon extending Reissner's 
and Griffith/Taylor's approaches, respectively, where the shear center location based 
upon Griffith/Taylor’s approach is clearly dependent upon (v) for thin sections, whereas 
the extended Reissner-based prediction is independent of (v). The circular symbol 
represents the closed form thin-plate prediction of Reissner (1989) for triangular isotropic 
cross-sections. These results illustrate: (1) for very low aspect ratio triangles the 
difference in the two approaches for locating the shear center can be profound, (2) the 
thin plate solutions of Reissner (1989) ar e valid for only a very small aspect ratio range 
(b/a< 0.04), (3) both shear center predictions converge as the aspect ratio approaches 
that of an equilateral triangle (b/a= 1.155), then (x s = x s *=0), and (4) for large aspect 
ratios, both shear center predictions are nearly coincident and they converge to the 
cross-section mid-length (a/2). This occurs because for triangular cross-sections with 
large aspect ratios, the out-of-plane flexural warping function approaches that of a thin 
rectangle cross-section. 

In a second study, the shear center location was determined for a slender 
cantilever beam (L/a= 10) with a thin triangular cross-section (b/a=.01) composed of off- 
axis unidirectional high-strength graphite/epoxy fibers (-90 0 <a<90 0 , J3=0°). In Fig. 9, 
both shear center locations are presented for the beam root and beam tip as a function of 
orientation angle (a). At the beam root, the extended Reissner based approach is 
independent of (a) whereas the extended Griffith/Taylor approach is highly dependent 
upon (a). This is expected, basedjjpon the above results for an isotropic triangular 
section which showed that the Griffith/Taylor solution is sensitive to material cross- 
coupling effects for thin sections. At the beam tip, both approaches for locating the shear 
center produce identical results when the orientation angle is either close to 0° (-10°<a< 
10°) or close to 90° (-70°<a<-1 1 0°, 70°<a<1 1 0°). Outside of this range, the two 
approaches produce tip shear center locations that can lie well outside of the cross- 
section shape, where the extended Reissner approach is much more conservative. 

in a third investigation, the sensitivity of the shear center location with varying (V6) 
was determined by studying a slender cantilever beam (L/a= 10) with a triangular cross- 
section (b/a=.1) composed of off-axis unidirectional high-strength graphite/epoxy fibers 
(o-O 0 , -90°</k90°). Since (V 4 =vs= 0 ) the shear center location is a cross-section 
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property that is independent of beam-length. In Fig. 10, both shear center locations are 
presented as a function of orientation angle ( p ). In the x-direction, the Griffith-Taylor 
based prediction is farther from the centroid and more sensitive to (/3) than the extended 
Reissner approach, whereas in the y-direction, the extended Reissner approach is 
clearly dependent upon (/3) and the Griffith-Taylor based prediction is virtually zero. 

NACA-001 2 Cross-Section 

A final study was performed to investigate the variation of the shear center location 
with material orientation angle (a) in typical composite general aviation aircraft wings 
and helicopter blades. These structures are approximated as homogeneous cantilever 
beams having a NACA-001 2 airfoil cross-section (Fig. 3.c), where a second coordinate 
system ( x,y ) is introduced with the origin at the leading edge. Two beam-lengths were 
considered; (L/c= 3) for typical general aviation aircraft wings and (Uc= 20) for long 
slender helicopter blades, where (c) is the cross-section chord. The airfoil cross-section 
is discretized using 95 points on the boundary (i.e. 95 triangular subregions), based 
upon the mathematical definition of (Abbott and Von Doenhoff, 1959), and the section 
centroid is located at (0.42067c). Each of the three local deformation functions are 
modeled using a 9 th order power-series polynomial (161 total unknown coefficients) and 
the numerical integration was performed exactly using a 52-point Gaussian Quadrature 
formula (Dunavant, 1985). 

In Fig 1 1 , both shear center locations are presented for the beam root and beam 
tip (L/c= 3,20) as a function of orientation angle (a). At the beam root the results are 
similar to the above triangular cross-section study (Fig. 9), where the extended Reissner 
based approach is nearly independent of {a), whereas the extended Griffith/Taylor 
approach is slightly more dependent upon (a). At the aircraft wing tip (Fig. 1 1, center 
region), it is observed that: (1) the two shear center definitions are in near perfect 
agreement when (-5°<a<10°) and (60°<a<1 10°), (2) the extended Reissner definition of 
the shear center can be ahead (x<0) and outside of the airfoil section if (-8°<a<-52°) or 
behind (x>c) and outside the airfoil section if (15°<a<42°), whereas the extended 
Griffith/Taylor definition of the shear center can be ahead (x<0) or behind (x>c) the airfoil 
section if (-8°<a<-60°) or (12°<a<55°), respectively, and (3) the maximum distance that 
the shear center can be located either ahead or behind the wing-tip section occurs for the 
extended Reissner definition at (ct=-30°, x*sL = '*5c) and (a=30°, x* s l = 1-2o), respectively, 
whereas the extended Griffith/Taylor definition has (a=- 38°, x s l=-1.1c) and (o=36o, 
x s l=1.7c), respectively. 
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For the tip section of the helicopter blade (Fig. i 1 , upper region), the shape of the 
curves represent an amplified version of the wing section. Thus, the shear center 
location is well outside of the tip cross-section for most values of (a), where the maximum 
distance that the shear center can be located either ahead or behind the wing-tip section 
occurs for the extended Reissner definition at (a=-30°, x* s l=-5c) and (a=30°, x* s l=6.5c), 
respectively, whereas the extended Griffith/Taylor definition has (a=-38°, x s l=-9c) and 
(o=36°, x s l=9c). respectively. 


Conclusions 

The behavior of a tip-loaded anisotropic cantilever beam with an arbitrary cross-section 
is studied using Saint-Venant's semi-inverse method along with a power series solution 
for the local in-plane and out-of-plane deformation warping functions. The power series 
coefficients are determined by solving a set of variationally derived linear algebraic 
equations. Using the resulting three-dimensional displacement solutions, the shear 
deformation associated with applied tip loads is investigated as well as the shear center 
location. Two different definitions of the shear center are presented for anisotropic 
beams by extending the classic approaches of Griffith/Taylor and that of Reissner. Both 
of the extended definitions reveal the linear dependency of the shear center location with 
beam-length, where the extended-Reissner prediction is much closer to the centroid then 
the extended G riff it h/Tay I o r prediction. Numerical results are presented for three different 
cross-sections and two different materials. For elliptical cross-sections, it was shown that 
the calculated power series coefficients were in exact agreement with existing elasticity 
solutions for anisotropic beams over a wide variety of cross-section aspect ratios. For the 
triangular cross-sections, it was shown that the calculated power series coefficients 
represent a "best-fit" to the infinite series of transcendental functions and the warping- 
related section properties (shear center, torsion constant, shear deformation) converge 
quickly with increasing power series order. Moreover, three studies were performed to 
illustrate the sensitivity of the shear center location with cross-section aspect ratio, 
material definition, fiber orientation, and beam-length. A final investigation was 
performed to study the length-dependency of the shear center in composite general 
aviation aircraft wings {L/c= 3) and helicopter blades ( 1 / 0 = 20 ). At the beam root, the 
extended Reissner approach is nearly independent of material orientation angle, 
whereas the extended Griffith/Taylor approach is dependent. At the aircraft wing tip, it is 
observed that the two shear center definitions are in near perfect agreement over a small 
range of orientation angles and the shear center can be located either ahead or behind 
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the wing-tip section. For the helicopter blade tip section, the shear center location is well 

outside of the tip cross-section for most values of orientation angle. 
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The matrix L J is defined as; 
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Integrating the matrix of (A1) over the beam length results in 



Table 1: 


Material properties for unidirectional high-strength graphite/epoxy fibers. 


£11 

145 

GPa 

£22 = £33 

10 

GPa 

Gi2 = Gl3 

4.8 

GPa 

vi 2 = vi 3 

0.250 


V23 

0.400 



2zn 


Table 2: Calculated section properties of an isotropic cantilever beam with a thin 

triangular cross-section as a function of power series order (b/a=0.1 , 
Ua= 10, v=0.333) 


Polynomial 

order 

matrix 

size 

xs/a 

xVa 

GJ (10-3) 

u shear OO' 3 ) 
Ototal 

^shear 00*5) 
Vtotal 

2 

14 

.2718 

.2923 

2.168 

4.640 

-24.26 

3 

26 

.1398 

.2035 

2.118 

5.815 

4.594 

4 

41 

.1441 

.2053 

2.060 

5.815 

5.049 

5 

59 

.1481 

.2070 

2.007 

5.815 

4.583 

6 

80 

.1515 

.2082 

1.969 

5.815 

4.671 

7 

104 

.1530 

.2089 

1.943 

5.815 

4.690 

8 

131 

.1536 

.2093 

1.938 

5.815 

4.700 

9 

161 

.1536 

.2093 

1.936 

5.815 

4.701 
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Table 3: Calculated section properties of a generally anisotropic cantilever beam with a thin triangular cross-section 

as a function of power series order (b/a=0.1 ,L/a=10, a=/3= 30°). 
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Fig. 2 Line of shear centers in an anisotropic beam. 






a (degrees) 


Fig. 5 Nondimensionalized in-plane coefficients (a01.a21.a03) as a function of 

orientation angle (a), (/J=0) and cross-section aspect ratio ( b/a ) for an elliptical 
cross-section ( — exact, 0 current approach). 
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Fig. 7 Nondimensionalized in-plane coefficients (c 01 .c 21 .c 03 ) as a function of 

orientation angle (a), ((3=0) and cross-section aspect ratio (b/a) for an elliptical 
cross-section ( — exact, 0 current approach). 





Fig. 8 Variation of the shear center location in an isotropic triangular cross-section as a 
function of ( b/a ) and (v), ( — Griffith/Taylor approach, — adaptation of 
Reissner approach, o Reissner prediction). 
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Fig. 9 Variation of the shear center locations at the beam root and beam tip of a thin 
composite triangular cross-section ( b/a ) as a function of orientation angle (a), 
{fi=0°). 





Fig. 10 Variation of the shear center locations at the beam root and beam tip of a thin 
composite triangular cross-section ( b/a ) as a function of orientation angle (/?), 
(0=90°). 
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Fig. 1 1 Variation of the shear center locations at the beam root and beam tip of a 
composite aircraft (L/c=3) and a composite helicopter blade (L/c= 20) as a 
function of orientation angle (ct), (p= 0°). 
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Chapter 4 Exact Cross-Section Warping Functions for Generally Anisotropic 
Beams Having Solid Elliptical Cross-Sections 

. . . ... , by 

C. A. le and J. B. Kosmatka 

Department of Applied Mechanics and Engineering Sciences 
University of California, San Diego 
La Jolla, CA 92093 

ABSTRACT 

The St Venant displacement distributions are developed, based upon the theory of 
elasticity, for a tip-loaded homogeneous cantilever beam having an elliptical cross- 
section and rectilinear anisotropy. These distributions are found by integrating the strain 
distributions, where the local in-plane deformation and out-of-plane warping of the cross- 
section are exactly determined. A definition for the 'anisotropic shear centef is presented 
by extending the classical definition for isotropic beams. The additional transverse beam 
displacement associated with shear deformation is determined for applied extension and 
flexure loads. Numerical results are presented which show (1.) the anisotropic shear 
center location is linearly dependent upon beam length and can be located outside the 
cross-section, (2.) shear deformation can actually be negative for certain beam aspect ra- 
tios and material definitions, and (3.) the local cross-section deformations and the trans- 
verse shear stress distributions bear no resemblance to their isotropic counterparts. 

INTRODUCTION 

The elastic stress and displacement distributions of isotropic cantilever beams 
subjected to tip loads (i.e., extension, bending, torsion, and flexure) has been exhaus- 
tively investigated by making use of Saint-Venant's principle in the formulation of the 
boundary-value problem. Closed-form displacement and stress solutions exist for simple 
(elliptical) cross-sections, series solutions exist for slightly more complex (rectangular, tri- 
angular) solutions, and approximate solutions based upon the application of the Ritz 
method exist for arbitrary cross-sections. Detailed examples of these solutions can be 
found in many texts covering the theory of elasticity (for example, (Sokolnikoff, 1956) or 
(Timoshenko and Goodier, 1970)). 
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Conversely, the study of generally anisotropic cantilever beams subjected to tip 
loads has received far less attention. In Lekhnitskii’s monograph (1963), the stress and 
displacement distributions were formulated in terms of known quantities (geometric and 
material properties) and unknown functions which represented the local in-plane defor- 
mation and out-of-plane warping of the cross-section. A solution procedure for determin- 
ing the stress distribution was presented based upon the use of Airy and Prandtl stress 
functions, where numerical examples include beams having an elliptical (closed-form), 
rectangular (series), or an arbitrary cross-section (approximate). But no results were pre- 
sented for the displacement distributions. Recently Kosmatka and Dong (1991) devel- 
oped an analytical (finite element based) model for determining the complete displace- 
ment and stress distributions of a homogeneous prismatic anisotropic beam with an arbi- 
trary cross-section, by solving the two-dimensional boundary problem in terms of the local 
in-plane deformation and out-of-plane warping of the cross-section. Numerical results 
used these calculated displacement and stress distributions to study the beam behavior, 
determine important section constants, and show that the shear center location is linearly 
dependent upon beam length. 

In the current paper, we will develop the complete St. Venant displacement and 
stress distributions, based upon the theory of elasticity, for a tip-loaded homogeneous 
cantilever beam having an elliptical cross-section and rectilinear anisotropy. The dis- 
placement distributions are found by integrating the strain distributions calculated by 
Lekhnitskii (1963), where the local in-plane deformation and out-of-plane warping of the 
cross-section are exactly determined. The resulting displacement solutions are used to 
develop a definition of the 'anisotropic shear centet which involves extending the original 
work of Griffith and Taylor (1917) to homogeneous prismatic anisotropic beams. 
Moreover, the additional transverse beam displacement associated with shear deforma- 
tion is determined for applied extension and flexure loads. It is also shown that the shear 
deformation is zero for applied bending and torsion of homogeneous anisotropic beams. 
Numerical results are presented to show how material anisotropy effects (1) the local in- 
plane deformation and out-of-plane warping, (2.) the shear center location, and (3.) the 
stress distribution. 

This model will be useful to developers of refined (or higher-order) one-dimen- 
sional theories for anisotropic beams, who wish to include the local in-plane deformation 
and out-of-plane warping of the cross-section as part of the displacement kinematic field 
In order to get exact three-dimensional stress distributions. 
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DISPLACEMENT DISTRIBUTIONS 


Consider a cantilevered beam of length (L) with an elliptical cross section of area 
(A) composed of an anisotropic homogeneous material (see Fig. 1 ). A Cartesian coordi- 
nate system ( x,y,z ) is defined on the beam where (x) and (y) axes are coincident with the 
principal (i.e., major and minor) axes of the root cross section and ( 2 ) is coincident with the 
line of centroids. Displacements fn the (x,y,z) directions are defined as (u), (v), and (w), 
respectively. At the root end (z=0), the beam is constrained by fixing the displacements at 
the centroid and the rotations about the centroid as follows 


u - v = w = 0 , 


du dv _dv 
9 z 9 z 9x 9y 


(la-/) 


At the beam's free end, tractions are applied that are equivalent to a general force (P) with 
flexural (P x , P y ) and extensions (P z ) components and a general moment ( M) with bend- 
ing (M*. My) and a torsional (A^) components. It is further assumed that the lateral surface 
of the beam is traction-free. 


The stresses and strains at a point can be written in array form as: 

{o) T = { Oxx, Oyy, Ozz , Tyz, T zx, *xy } , 

{e} 7 ’= { £xx, £yy, £zz, tyz, tyx, tyy ) . 


where the strains are related to the displacement components by: 


du 

£xx = — , 

9x 

9v 9w 
Yyz = — + — , 

9 z 9y 


9v 

eyy = Yy 

9 W du 

Jzx — -r— + — 

9x 9 Z 


£zz = 


9w 


9z 

du 


dv 


Yxy - tt— + ■_ 

9y 9x 


( 2 . 0 / 7 ) 


and the stress and strain components are related by the constitutive relations of a lineariy- 
anisotropic hyperelastic material: 


(£) = [S]{o). 


( 2 ./) 


Here [S] is a fully populated symmetric matrix with 21 distinct coefficients Sjj (i,j- 1-6). 


The behavior of the beam will be studied as two independent cases. The first case, 
which involves generalized plane strain behavior of the beam, is associated with the ex- 
tensional force (P z ), the bending moments (Mx, My), and the torsional moment (M z ). The 
second case is associated with the applied flexural forces (P x , Py), where it is assumed 
that the stress (czz) varies linearly in the z-direction and the remaining five stress compo- 
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nents can be nonzero. 


Case I: Generalized Plane Strain 


For a homogeneous anisotropic cantilever beam with an elliptical cross-section 
subjected to an applied extensional force (P z ), bending moments (M x , My), and/or tor- 
sional moment ( Mz ), the stress distributions are given (from Lekhnitskii (1963)) as: 


Oxx = ®yy — Txy — 0 , Ozz — ^ + j 


y My 


XX 


i 


Xyz = 


Mr 


yy 


2/, 


Tzx 


yy 


m 2 

21 


(3.a-t) 


XX 


where (a) and ( b ) are one-half of the major and the minor dimensions of the cross-section, 
respectively, (see Fig. 1), A is the cross sectional area given as (itab), and l xx and lyy are 

the area moments of inertia about the (x) and (y) axes given as (na&IA) and (^bM), re- 
spectively. These stress components, which are independent of ( z ), are equal to those of 
an identical isotropic beam. 


The displacement components ( u,v,w ) are determined by substituting Eqns. (3 .a-f) 
into (2./) and integrating using the standard technique (see Sokolnikoff (1946) or 
Timoshenko (1970)), where the constants of integration are determined using the bound- 
ary conditions (1 .a-f) at the root. Thus, 



+ (^-z 2 ) - jj £- ( 2 S 15 X +S 56 .K +Sssz)y +jj £ -{S-\4X 2 -S24y 2 -S44yz) . 

d.clyy\d I 4 l xx **iyy 


v * v ?y} + 7^j-( v s x:! + 2v 2xy +«>«)- 2 U a 

" “ ( 4 .b) 

~^r[% z2 V 77 t (Sisx2-S 2s y2 + S 55 «) + + 2S 24 y +S442)/ , 

2E l xx \2 I A I XX 4lyy 




Ml 

Ixx 


. , (S55* + §45/ + 2$35Z)y + -jjH S45X + Suy + 2§34Z)x , 

4 Ixx 4 lyy 


(4.c) 


in which E (=1/533 ) ' s Young's modulus in the z-direction and v/ are the cross-coupling 
coefficients defined as vj=-Sj3/S33, where v j and v 2 are the usual Poisson coefficients, 
and v 4 , V5, and v 6 express the three-dimensional extension-shear coupling in the 
anisotropic beam. The above displacement results are in agreement with the plane-strain 

4.4 


solutions developed by Lekhnitskii (1963) for an anisotropic beam. One can see that the 
application of an extension force (P z ) will produce beam extension, shear deformation (as 
a result of v 4 and V5), and deformations within the cross-section plane composed of 
Poisson contractions (from v 1t v 2 ) and shear (v 6 ). Applying a bending moment {M x , My) 
will produce beam bending as well as beam twisting (as a result of v 4 and v 5 ), in-plane 
cross-section deformations which include the formation of an anticlastic surface (v-j, v 2 ) 
and shearing (v 6 ), and out-of-plane cross-section deformation that resembles torsion-type 
warping (bi-linear function) as a result of anisotropic material coupling (v 4 ,V5). Finally, 
applying a torsion moment ( M z ), will produce beam twisting as well as beam bending (as 
a result of v 4 and V5) and shear deformation, in-plane cross-section deformations which 
include shearing and the formation of an anticiastic-type surface, and out-of-plane cross- 
section deformation associated with torsion-warping (bi-linear function). 

Case II: Applied Flexure Load 

The behavior of the homogeneous anisotropic cantilever beam subjected to ap- 
plied flexural load (P y ) can be studied using the following stress functions (from 
Lekhnitskii (1963)) 



where the coefficients (£?/, hi, 2,3 ), which are functions of the cross-section aspect ratio 
and the material definition, are determined by solving a set of linear algebraic equations. 
These three equations are presented in the Appendix (Eqns. A1-14). The stress compo- 
nents written in terms of the stress functions are 


Oxx = 


CTzz = - 


-.2 

3 Y*i 

dy 2 

P y 

2 1 XX 


Cyy = 


-.2 

a v*i 


3x2 

2 y[L - z) + \6 xy + v 4 y 2 + 




1 1X 


+ VI Gxx + VfcOyy + V4T yz + V$T zx + MS'Txy » 

(6 .a-f) 

3yg = 3 v^i 

dy y dx dy 


By substituting the stress functions (Eqns. 5.a,b) into (6 .a-f), one can easily sees that: (1) 
the in-plane stresses (<t xx , c yy , r xy ) are nonzero as a result of anisotropy and only func- 
tions of the cross-section coordinates, (2) the normal stress (<r zz ) contains the classical 
strength of material term ( -Py ( L—z)y/l xx ) an ^ remaining terms of (oiz) ar ® inde- 
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pendent of (z) and produce a zero resultant when integrated over the cross-section area, 
and (3) the shear stresses (r xz , Ty Z ) contain the classical isotropic terms and additional 
terms associated with anisotropy which are both independent of (z). 


The displacement components (u,v,w) are determined by first substituting Eqns. 
(5 a-/) into (6.a-f), follow by substituting the resulting stresses into (2./) and then integrating 
the strains, where the constants of integration are determined using the boundary condi- 
tions (1 .a-/) at the root. Thus, 


u = 


2 El 


XX 


(2 vi x + v 6 y) y{z- L) + ^y^z - 2L)j 


+ Pyjaio x+aoi y+a3o x 3 + a 2 i x 2 y+ai 2 xy 2 + ao3y 3 J , 


v' - - 2^if ( z ‘ 3L )' f xz ( z ■ 2L ) + ( vz y 2 ' Vl * 2 ) ( z ' L )) 

(7.a-c) 

+ P y |b 10 x+boi y+h ox 3 + b 2 i x 2 y + bi 2 xy 2 + bo3 y 3 | . 


w = -j^{{vsx+v A y)y{z-L)-yz{z-2L))-^yz 

+ Pyle 10 x+ CDI y + C 30 X 3 + CSi X 2 y+ C 12 xy 2 + PD3 y 3 ) . 


where the coefficients ay, bij, and a\ (i,j = 0-3) are presented in the Appendix (Eqns. A.15- 
32) and the subscripts (/) and (j) refer to the order of the polynomial (x^). The above 
displacement distribution exactly describes the extension, bending, and twisting of the 
anisotropic cantilever beam as a result of an applied flexure force (P x , P y ), as well as local 
in-plane cross-section deformations (all terms associated with ay and £*j) and out-of-plane 
cross-section warping (all terms associated with cyj. These solutions are in exact agree- 
ment with the results of (Kosmatka and Dong (1991)), where their model uses the finite 
element method to approximately determine the local cross-section deformations for an 
arbitrary cross-section. 


ANISOTROPIC SHEAR CENTER 


The ’ shear center ' for a prismatic beam composed of an isotropic material is a 
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property of the cross-section and independent of beam-length. It is located using the 
classical definition (Griffith and Taylor (1917)), 'the load point that produces a zero mean 
value cross-section twist (i.e., zero twist about the centroidal axis)'. In fact for doubly- 
symmetric cross-sections (including the current elliptical cross-section), the shear center 
of an isotropic prismatic beam is located at the centroid. For a prismatic beam composed 
of a generally anisotropic material (V4, V5*0), the classic definition is not applicable and 

we will develop a more appropriate definition. 

The twist about the centroidal axis (x=y=0) for the cantilever beam subjected to a y- 
direction tip flexure load (P 2 ) offset from the centroid by an amount (x s ) is calculated using 

0 = (8 a) 

2\ax By) 

where (u, v) are found by combining Eqns. (4 .a,b) and (7 .a,b) with (P y =P 2 ) and (A4 2 =P 2 x s ); 

6 = ^2_(S35(2L - z) + x s (s 44 (ff + S55)}* ■ M 

The single-underlined term, which is associated with anisotropic *bend-twist' coupling 
(Kosmatka and Dong (1991)), is a quadratic variation of the cross-section twist as a result 
of the applied flexure force (P y =P 2 ), whereas the double*underlined term represents a 
linear variation as a result of the applied torsion moment (A/f z =P 2 x s ). Since the cross-sec- 
tion twist varies quadratically along the beam-length, it is not possible to achieve zero 
twist along the entire beam-length. Instead the best that one can achieve is zero twist at 
two locations: the beam root (z=0) and the x-y plane that contains the applied load. Thus 
the x-direction 'anisotropic shear center ' location for the free-end of the cantilever beam 
is found by substituting (z=L) into Eq. (8 .b) and finding the value of (x s ) that produces a 

zero cross-section twist (6=0): 



where the location depends upon the beam-length and reduces to the classic definition 
(centroid) for isotropic materials. The variation of the cross-section twist over the beam- 
length is found by substituting Eq. (9) into (8.b) 

(10) 

where the twist is zero at the beam root and tip (load plane) and reaches a maximum at 
the beam mid-length (Fig. 2). 


4.7 


The y-direction component of the shear center can be found in a similar fashion by 
applying an x-direction flexure force (P^) that is offset form the centroid by an amount (y s ). 
Finally we formally define the 'anisotropic shear center ' as; "the load point that produces 
general bending and twisting of an anisotropic beam with zero mean value twist (zero 
twist about the centroid) in the cross-section plane that contains the applied load" 

SHEAR DEFORMATION 


An examination of the displacement components of Eqns. (4.a-c) and (7 .a-c) re- 
veals that the deformation of the centroidal axis (x=y=0) agrees with the standard strength 
of materials solutions, but that the additional displacement associated with shear defor- 
mation is not included. This occurs because our original root boundary conditions (Eqns 
(1 .a-f)) assume that the slope of the centroidal axis is zero ( du/dz= BvlBz- 0) as opposed 
to fixing the rotation of the root cross-section ( dwldx = dw/dy = 0). This additional dis- 
placement associated with shear deformation can be included by simply rotating the de- 
formed beam so that slope of the deformed root cross-section at (x=y=0) is coincident with 
the x-y plane, and thus the deformed centroidal axis will have a nonzero root slope (see 
Timoshenko and Goodier (1970) for further details). These rotation angles are equal to 
the shear strains (ft z , fy z ) at the beam root centroid (x=y=z=0). Calculating the rotational 
angles using Eqns. (4.a-c) and (7.a-c): 


= Yyz 
x-y=z = 0 


■ P >< C01 >' P2 (^I ’ 


4>y ~ Yxz 

x-y=z = 0 


- P^ C , 0 )-P^J , (11.3,6) 


where shear deformation for a homogeneous anisotropic beam occurs as a result of ap- 
plied extension and flexure, but not applied bending or torsion. The final form of the dis- 
placements distributions that includes shear deformation is found by modifying the dis- 
placements (u,v,w) of Eqns. (4.a-c, 7.a-c) in the following manner: 

Z? = u + z <py , ~v = v+zfr , 'w = w-yfo- xty. (12.a-c) 


NUMERICAL EXAMPLES 


In this section, three detailed numerical examples are presented to show how ma- 
terial orthotropy or anisotropy effects the local in-plane deformation and warping of the 
cross-section, the shear center location, and the shear stress distribution. Only the behav- 
ior of the beam as a result of an applied flexure force (P y ) will be studied, since beam be- 
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havior as a result of extension (P z ), torsion (M z ), and bending {Mx, My) has been studied 
in detail by Lekhnitskii (1963) and Kosmatka and Dong (1991). 

The beam is assumed to be composed of a single set of unidirectional high 
strength graphite/epoxy fibers, which is a transversely isotropic material that has five in- 
dependent constants (see Table 1) with a sixth constant that is equal to G 23 - 
E22/(2(1+V23)). To achieve orthotropic and/or anisotropic beam behavior, the material ref- 
erence frame (1,2,3) is oriented relative to the beam Cartesian coordinate frame (x,y,z) 
using reference angles (0) and (a), which are defined as rotations about the positive z- 
axis and the positive material 3-axis, respectively. See Fig. 3, where the transformation 
relation between (1 ,2,3) and (x,y,z) is equal to 


M 

2 


sin(a) cos(P) 
cos(a) cos(P) 

sin(a) sin(P) 
cos(a) sin(P) 

cos(a) 
- sin(a) 

( : ) 

i 3 j 


- sin(P) 

cos(p) 

0 



The resulting 21 unique material compliance coefficients (Sy, i,j=1-6) are determined us- 
ing standard transformation techniques (see Lekhnitskii, 1963). 

E xam pte_l 

In this example, the local cross-section deformation and the shear center location 
are studied as a function of material orthotropy and cross-section aspect ratio. Material 
orthotropy is introduced by rotating the fiber set in the x-z plane using the reference angle 
(a), while holding (/3=0). Note that, in this situation, Sj4=Sj6=0 (^1-3,5) and V4=ve=0. 


The coefficients (£?/ , i = 1 ,2,3) are determined using Eq. (A1) as B-\ - S 3 = 0 and 



Sjj 5 - 4S13S33 + 2 S 33 S 44^) 2 
Sfj 5 - S 33 S 55 - 3S33S44 ^) 2 


(14) 


The local deformation within the cross-section (x-y) plane and warping out of the 
cross-section plane as a result of an applied flexural force (Py) is characterized in the u, v, 
and w directions, using Eqns. (7.a-c), as simply 

P^aoiy+azix^+aosy 3 ), P>|biox+bi2xy 2 +b3ox 3 j, Pj|ooiy + P2i x2 / + ®03y 3 )* (15.a-c) 

where the other deformation coefficients (ay, by, cy) are zero for ((}= 0). In Figs. 4 through 
1 1 the nonzero coefficients are presented as a function of orientation angle (a) and 
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aspect ratio ( b/a ), where they have been nondimensionalized using either (E11A) or 
(Ei i /xx)- For this material, all of the deformation coefficients are independent of aspect 
ratio (b/a) for a=8° and 60°. The coefficients aoi and bio. which define the uniform in- 
plane cross-section shear, are presented in Fig. 4, where positive values can be obtained 
for thin cross-sections (b/a=0.1) with 18°<cr <39°. In Figs. 5 and 6, the coefficients (a2i 
and 630) which are associated with (x 2 ) variations in (fry) are presented. These two coef- 
ficients, which have relatively large magnitudes and opposite signs, represent a signifi- 
cant amount of in-plane cross-section deformation associated with relatively low (ft y ). 
The remaining coefficients (ao3 and bi2), which are associated with (y 2 ) variations in 
(fty), are presented in Figs. 7 and 8. Although these coefficients are generally smaller 
than a2i and £>30, they are associated with a larger component of (ft y ). For thin cross- 
sections (b/a=0.1), ao3 reaches a positive maximum at 60°, bi2 reaches a negative maxi- 
mum at approximately 33°, and the resulting component associated with (y 2 ) variations in 
(tty) reaches a positive maximum at 60° and a negative maximum at 30°. Whereas for 
thicker cross-sections (b/a=1.0,10.0), the component associated with (y 2 ) variations in 
(tty) is always positive and reaches a maximum at 45°. 

The coefficient (coi), which is proportional to the amount of y-direction shear de- 
formation associated with P y (see Eq. 11. a), is presented in Fig. 9, where for thin cross- 
sections (b/a=0.1) this term is becomes negative for 18°<a <39°. The out of plane warp- 
ing constant (C21), which is always negative for thin cross-sections (b/a=0.10), is pre- 
sented in Fig. 10. The remaining out of plane warping constant (C03) is presented in Fig. 
1 1 , where it is always negative for all but the thinnest cross-sections. 

The 'anisotropic shear center location ' of the beam tip cross-section (z=L) as a re- 
sult of an applied transverse tip load (P2) can be determined by substituting the definition 
of the compliance coefficients with (/MD) into Eq. (9) 

x /I (sT &7H n(2a) * fe' a^ 1 + 2V|3> ~ 

s 2 cos 2 (2«) | (1 ( Eh , 2 v . |3) s/n 2 (2«) t foV2 IcosHo) ( 2 (1+vs 3 )s/n 2 (g) j 
\ G13 E22 E11 a ' G13 E22 I 

where the results are presented in Fig. 12 as a function of orientation angle (a), aspect ra- 
tio (b/a), and normalized to the beam length (L). For (a=0°) and (80°<a<90°) the shear 
center is located at the centroid (i.e. isotropic shear center). For (0°<a<80°) P the 
anisotropic shear center is offset from the centroid, with the maximum offset occurring 
near (a=40°) for thin cross-sections (b/a=0.1) and approaches (30°) for circular (b/a= 1) 
and thicker cross-sections (b/a=5). As an illustration, the anisotropic shear center of a 
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slender beam (L/a= 10) with a thin cross-section (b/a=0.10) and (a=40°) can be found us- 
ing Eq. (16) or Fig. 12 as x s = 4.8a, which is well outside the cross-section. Furthermore, 
doubling the beam length (L/a= 20), will double the offset to the shear center location 
(x s =9.6a). A simplified form of lie shear ce nter location for small orientation angles (a) 
can be found by taking a Taylor series expansion of Eq. (16) 



(17) 


Example 2 

A second example is presented to explain why the added displacement associated 
with shear deformation (cq-\L) can be negative for a thin (b/a= 0.1) tip loaded (P y ) can- 
tilever beam composed of off-angle unidirectional graphite/epoxy with 18°<a <39° and 
/3=0°. See Fig. 9. Thus the inclusion of shear deformation makes the beam stiffer, as op- 
posed to more flexible as one normally sees in plate bending. From Eqn. (A.28), 

Coi = Aw(^ + 32)» (1®) 

where this unusual situation occurs whenever ( Bz < -2/A), since $44 which equals S44 for 
(cej*0, /3=0) is always greater than zero. Substituting Eqns (5 .a,b) into (6.d), the transverse 
shear stress (T yz ) can be written as 

where the shear stress (r yz ) at the centroid (x=y=0) will be negative whenever B 2 < -2/A. 
The transverse shear stress (r yz ) distribution for two slices of a beam element (dz) is pre- 
sented in Fig. 13, where the stresses are small and negative near the centroid and large 
and positive near the outer edges (x=±a) so that the integral of (t yz ) over the cross-section 
will equal (P y ). This situation is typical for elliptical cross-sections and will not occur for 
rectangular (plate-like) cross-sections, since (T yz =0) on the upper and lower surfaces of 
the rectangular (plate-like) cross-section to satisfy the traction free boundary conditions. 


An expression can be developed to define the bounds on this unusual condition by 
substituting Eq. (14) into (20) and rearranging 

Itif < 3S| 5 - 2S33(S55 + 2 Si a) ( 20 ) 

l a J 4S33S44 




For an isotropic beam, the above expression reduces to (b/a) 2 < -1/(2(1+v)) and thus this 
condition will not occur. For an anisotropic beam composed of off-angle uni-directional S- 
glass/epoxy (Eii - 55 GPa, E22=P33=16 GPa, Gi2=£»i3=7.6 GPa, v-j 2 = vi 3=0.28), the 
normalized shear deformation coefficient (qqi) is presented in Fig. 14 for a broad range of 
(a ) and (b/a). It can be seen that this unusual condition does not occur and thus the in- 
clusion of shear deformation will only make the beam more flexible. This can be traced to 
the fact that S-glass/epoxy is closer to an Isotropic material than the graphite/epoxy of ex- 
ample 1. It is interesting to note that for S-glass/epoxy, all of the deformation coefficients 
are independent of aspect ratio (b/a) for a =22° and 55° and the width between these two 
points is smaller than that of the graphite/epoxy. Finally, as one chooses a material that 
closely resembles an isotropic material (G 12 * 5ii/(2(1+vi2)) ). then these two points dis- 
appear so that all curves are independent of (a). 

Example 3 

In this final example, the local cross-section deformation and the shear stress dis- 
tribution are studied for the case of a general anisotropic (a,/J *0) beam composed of 
graphite/epoxy (Table 1) and subjected to a flexural force (Py). For this general case, all 
three coefficients (8j, *=1-3) and all eighteen deformation coefficients (ay, by, Qj, i,j= 1-3) are 
nonzero. The local in-plane deformations and the out-of-plane warping of the cross-sec- 
tion are presented in Fig. 15.a,b for the case of (b/a= 0.5) with orientation angles of 
(a=30°) and Q3=60°). The local in-plane deformation (Fig. 15.a) is represented by a non- 
symmetric shape that is dependent upon the material orientation angles (a,/3). The out-of- 
plane warping (Fig. 15.b), which does not resemble the symmetric cubic-like warping of 
an isotropic beam (Sokolnikoff, 1956), is divided into six alternating regions of positive 
(solid lines) and negative (dashed lines) deformations. The effects associated with the 
shear deformation coefficients (c-io, coi) are not included in Fig. (15.b) because they pro- 
duce only a rigid cross-section rotation and no cross-section deformation. 

A quantitative vectorial scaled plot of the transverse shear stresses (f^x* Tyz ) for 
(a/b= 0.5) with orientation angles of (a=30°) and (0=30°) is presented in Fig. 16. From this 
figure, it is interesting to see that the stress distribution does not resemble that of an 
isotropic beam (i.e. a paraboloid with the maximum occurring at the centroid). Instead, the 
distribution is nonsymmetric with the maximum occurring on the outer edge. 

CONCLUSION 

The complete St Venant elastic displacement and stress distributions are devel- 
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oped, based upon the theory of elasticity, for a tip-loaded homogeneous cantilever beam 
having an elliptical cross-section and rectilinear anisotropy. The displacement distribu- 
tions are found by integrating the strain distributions, where the local in-plane deformation 
and out-of-plane warping of the cross-section are exactly determined. A definition for the 
• anisotropic shear center 1 is presented based upon extending the classical definition for 
isotropic beams. The additional transverse beam displacement associated with shear 
deformation is determined for applied extension and flexure loads. Numerical results are 
presented which show for the flexural behavior of an orthotropic (c& 0,f$=0) beam that the 
local in-plane deformation and out-of-plane warping are highly dependent upon fiber ori- 
entation and cross-section aspect ratio, where for two orientations the local deformations 
are completely independent of cross-section aspect ratio. The anisotropic shear center 
location is linearly dependent upon beam length and can be located outside the cross- 
section depending upon the fiber orientation. Moreover, the added transverse displace- 
ment associated with shear deformation can actually be negative for certain beam aspect 
ratios and material definitions so that the inclusion of shear deformation may actually 
make the beam stiffer, as opposed to more flexible. Finally, the local cross-section defor- 
mations and the transverse shear stress distributions bear no resemblance to their 
Isotropic counterparts. 
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Appendix 


The constants (S/, / ■ 1 , 2 , 3 ) are determined by solving 


'G11 

G12 G13 


Sl 

Hi 

G21 

G22 G23 


B 2 

- Hz 

. G31 

G32 G33. 


1 03 1 

1 h 3 1 


(A1) 


where the coefficents of [G] and [HJ are defined as: 

Gn = (30i 1 (^) 2 +3 022(aF + 2Ai2+A>6) . G 12 * ■ ( + + ^ ’ 

G 13 = ( 3yJts + (fe + M|f ) ■ (A<+3ft4af + & 6 > • 

G 22 = ( 30«(§f + fts ) , ^3 * • ( ^^(a) 2 ^ ’ 

GB1 = ^ 0325 + A*6 + ) > G32 = ' 2045 . 

G 33 = ( Aw^f + 3055 ) . 

(Si4 * 20u) - ^($56 - 056) - A24(aF} * 

1(S55 - £ 55 ) * 2 S 13 + Aw(aFj ’ 

1(S45 - 3^45) * $36 1 . 

S31 S$J 


(A 2 -A 10 ) 


and 


* = .Z 


(All- 13 ) 


Ay- % 


S33 


(A. 14 ) 


The x-direction deformation coefficients associated with an applied flexure load: 


aio = - 4 Si 


^-• 2S te + fHxHf! + ^ + 77?' 


^ = 4 S L |All + Ai2 


921 = • Si 


a 2 |3p 2 a 2 
4^6 




(A. 15 - 18 ) 


|a 2 b 2 


+ $> 


Ais 


Htl-S • 


4.14 


irr fp $' in if f" w 


(A 19-20) 



The /-direction deformation coefficients associated with an applied flexure load: 




The out-of-plane deformation coefficients associated with an applied flexure load: 

Cb, - • 4S,{^ + ♦ bJjU, | - feffe ) * ^ ^ . (A27-29) 



j Pis | (ks\ 

W a2 / 






Table 1 : Typical Material Properties of High-Strength Graphite/Epoxy. 


Eu 

145.0 GPa 

E 22 = e 33 

10.0 GPa 

Gi2= G|3 

4.8 GPa 

vi2 = v 13 

0.25 

v 23 

0.40 


4lfe 
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Fig. 1 A cantilever tip-! 






Nondimensionalized in-plane coefficients (aoi, 610 ) with (a), 03=0). 





a (degrees) 

Fig. 5 Nondimensionalized in-plane coefficient (a 2 1 ) with (a), (£=0). 






V* 




Nondimensionalized in-plane coefficient (ao3) with (a), (/J=0) 
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Fig. 9 Nondimensionalized shear deformation coefficient (cqi) with (a), (fi= 0). 










Fig. 1 1 Nondimensionalized out-of-plane coefficient (C03) with (a), 03=0). 





Fig. 12 Variation of normalized shear center location (x s ) with (a), 03=0). 
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Fig. 13 Transverse shear stress distribution (tyz) on two slices of a thin cross-section 
beam element. 
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Fig. 14 Nondimensionalized shear deformation coefficient (cqi) with (a), (/3= 0) for S- 
glass/epoxy. 


2 a 



Fig. 15 


a.) In-piane deformation, and b.) contour of out-of-plane warping of cross 
section with (£>/a= 0.5, a =30°, /J=60°). 
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Chapter 5 Transverse Vibrations of Shear-Deformable Beams 
Using a General Higher-Order Theory 

by 

J. B. Kosmatka 

Department of Applied Mechanics and Engineering Science 
• • ■ University of California 
San Diego, California 92093 


1. INTRODUCTION 

The fundamental vibration behavior of long slender cylindrical or prismatic beams 
can be studied using the classical Bernoulli-Euler beam theory. Attempting to use this 
theory to study either short beams or higher vibrational modes can lead to a significant 
over-prediction of the natural frequencies since both transverse shear deformation and 
rotatory inertia effects have been ignored. Timoshenko [1,2] developed a theory which 
allows one to study the vibrational behavior of shorter beams (and/or higher vibrational 
modes) by approximately accounting for both transverse shear deformation and rotatory 
inertia, where the cross-section remains undistorted and the in-plane stresses are zero. 
The resulting model is characterized by two differential equations of motion 
encompassing two independent variables; the transverse deflection of the neutral axis, v 
and the rotation of the cross-section measured about the neutral axis, 6. In addition, the 
model requires the determination of a well-known shear-correction factor (k), which is 
defined as the ratio of the average shear strain within the cross-section to the shear strain 
at the section centroid. Other researchers [3-5] have attempted to show that the 
magnitude of ( k ) should be adjusted for studying the higher mode vibrational behavior of 
beams because the dynamic shear strain distribution may differ significantly from the 
parabolic form of the static shear strain distribution. Cowper [6] presented a conceptual 
modification to Timoshenko equations by assuming that the transverse displacement, v, 
and the rotation, 6, represent the average cross-section values instead of the point-wise 
values. Thus, the ( k ) values could be determined, based upon static three-dimensional 
elasticity theory along with Saint-Venant's static flexure warping function, for a wide 
variety of cross-section shapes and Poisson ratios. A further study by Cowper [7] has 
shown, by comparison to plane-stress elasticity solutions, that the (/c) does not have to be 
adjusted for higher vibrational modes. 
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Goodier [8] was one of the first researchers to improve upon Timoshenko's beam 
theory by incorporating second-order stress effects associated with shear deformation 
and cross-section warping. While this model is limited to static behavior, it does provide 
valuable insight into beam behavior. Murthy [9] took a fresh look at the model of [1 ,2] by 
incorporating additional displacements to insure a zero in-plane stress state without the 
need for a correction factor ( k ). Stephen and Levinson [1 0] developed a second order 
model for studying beam vibrations that combines the best features of [2], [6], and [8] with 
two governing differential equations having two independent coefficients. The first 
coefficient is similar to Timoshenko’s (or Cowper's) shear correction factor that depends 
upon Saint-Venant's static flexure warping function, while the second coefficient 
depends upon the transverse shear stress state. 

Higher-order (or refined) shear-deformable theories [11-14] have been developed 
for beams with thin rectangular cross-sections that correctly account for the stress-free 
boundary conditions on the upper and lower surfaces as well as the parabolic shear- 
strain distribution through the thickness without the need for a shear correction factor ( k ). 
This is accomplished by expanding the axial displacement to include a cubic distribution 
through the thickness. This additional displacement is identical to the Saint-Venant's 
static flexure warping function for thin rectangular cross-sections, but for general cross- 
section shapes the correct expansion should be an infinite series of transcendental 
functions. These theories further assume that the in-plane cross-section stresses are 
negligible and that the cross-section does not deform in it's own plane. Thus, the 
possibility of dynamic deformation within the cross-section plane that occurs as a result 
of Poisson coupling with the out-of-plane cross-section stress distribution (anticlastic- 
type surface) is not included. The wave (or vibrational) behavior of a general prismatic 
bar with an arbitrary cross-section was studied by Aalami [15] using a Rayleigh-Ritz 
energy approach to the general three-dimensional problem. His numerical and 
graphical results clearly illustrate the presence and importance of both the out-of-plane 
shear deformation and in-plane (anticlastic) deformation for extremely short wave (high 
vibrational) modes (2 b/L = 3). Recently, le and Kosmatka [16] developed a static theory 
for a general cylindrical or prismatic beam that incorporates both out-of-plane shear and 
in-plane (anticlastic-type) deformation functions, where these functions are assumed 
known. In actuality, these functions can be determined exactly for simple cross-sections 
(rectangle, ellipse) by solving Saint-Venant's bending and flexure problems and 
approximately for an arbitrary cross-section by applying either a two-dimensional finite 
element approach [17] or a power series approach [18]. Numerical results show that the 
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calculated displacements and stresses are indistinguishable from elasticity solutions for 
a wide variety of beam loadings, boundary conditions, and cross-section shapes. 

The purpose of the current study is to develop a general higher order theory to 
study the static and vibrational behavior of beam structures having an arbitrary cross- 
section that utilizes both out-of-plane shear-dependent warping and in-plane 
(anticlastic) deformations. The equations of motion are derived via Hamilton's principle, 
where the full three-dimensional constitutive relations are used. In addition, a simplified 
version of the general higher-order theory is also presented for beams having an 
arbitrary cross-section that includes out-of-plane shear-deformation, but assumes that 
stresses within the cross-section and in-plane deformations are negligible. This 
simplified model, which is accurate for long to moderately short wave-lengths, offers 
substantial improvements over existing higher-order theories [11-14] that are limited to 
beams with thin rectangular cross-sections. Furthermore, the current approach will be 
very useful in the study of thin-wall closed-cell beams (for example, airfoil-type 
sections), where the magnitude of shear related cross-section warping is significant. 

A series of numerical results are presented that fully validate the current model 
with existing one-dimensional models as well as with appropriate elasticity solutions 
where available. The vibrational behavior of a simply-supported beam having either a 
rectangular or an elliptical cross-section is studied using the Ritz method for a wide 
variety of cross-section aspect ratios and beam (wave) lengths. Moreover, this problem 
contains a second (higher) spectrum of shear-dominant frequencies [19-22] and thus, 
the interaction of the current theory having dynamic in-plane deformations with the 
second frequency spectrum will also be investigated. 


2. GENERAL HIGHER-ORDER THEORY 

Consider a cylindrical or prismatic isotropic beam, of length L, having a general 
homogeneous cross section of area A. See Fig. 1. A Cartesian coordinate system 
( x,y,z ) is defined on the beam where x and y are coincident with the principal axes of the 
beam root cross-section and z is coincident with the centroidal axis. To eliminate 
complications associated with torsional vibrations, it is assumed that the centroidal axis 
and the elastic axis are coincident. The current development is further restricted to the 
study of vibrational behavior in the y-z plane only and the displacement relations are 
defined as; 


5.3 


(I.a-c) 


u(*,y,z,t) = M{z,t) Uo(x.y) , 
v(x,y,z,f) = v(z,f) + M(z,t) V’o(x.y) , 
w(x,y,z,t) = y 8{z,t) + Q(z,f) W^x.y) , 

where the components are defined in three groups. The first group contains the 
classical (or first-order) terms (v, yd), where v(z,f) represents the time-dependent y- 
direction displacement and 6{z,f) the time-dependent rotation about x axis. The second 
group includes the out-of-plane warping of the cross-section that is linearly proportional 
to the local time-dependent shear resultant (Q). The third group represents the in-plane 
deformation of the cross-section that forms an anticlastic surface and is linearly 
proportional to the local time-dependent moment resultant {M). The three functions 
(Uo,V 0 ,W 0 ), which are assumed to be known, can be determined exactly for simple 
cross-sections (rectangle, ellipse) by solving Saint Venant's bending and flexure 
problems [23], or approximately for arbitrary cross-sections using either a two- 
dimensional finite element (Ritz) approach [17] or a power series approach [18]. 

The time-dependent strain displacement relationships of the beam are defined as 


Bxx - U.x » Yyz - v ,z + w ,y » 

Byy = v,y , y xz = W' X + u t7 , 

£zz = W,z , Yxy = U,y +V ? .x , 


(2 .art) 


where exx, £yy, and ezz are the normal strains in the x, y, and z directions, respectively, 
and yyz, yxz, and yxy are the shear strains. The notation ( ) >x and ( ) % y refer to partial 
derivatives with respect to the x and y coordinates respectively. The general 
constitutive relations are given as 

(al = [C](e), 

where the stress and strain arrays are defined as 

(c}^"= { Cxx. Oyy, <Jzz, Tyz, *xz, Txy j . 

( Bxx, £yy,£zz> Yyz • Yxz, Yxy } , 

and the nonzero coefficients of the material stiffness matrix are given as 
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Cii = C 22 =Qj3 = (A + 2G) , 

C12 = C13 =C23 = A, ( 3 .d) 

C 44 = C 55 =C$q — G , 

Here, A and G are defined as Lame's constants. 

At this point, we assume that (Q(z,/), Z * 0). so that the displacement field of Eqns (l.a- 
c) can be expressed in terms of the kinematic variables ( v,6 ) only, instead of with a 
mixed form ( v,6,Q,M ). Clearly, for extremely short wave-length (higher vibrational) 
modes the current assumption that the warping deformation is only dependent upon the 
local bending moment and shear resultants, may lead to minor inaccuracies. For these 
types of modes, where the characteristic length is much less than the cross-section 
dimensions, one could further assume that the warping deformation is proportional not 
only to the local bending moment and shear resultants but also to their (first, second, 
etc.) derivatives as well. 

The local bending moment can be expressed in terms of the kinematic variables by 
making use of the cross-section equilibrium equation with 

= J yczzdA = El xx 0{z,t),z . ( 4 a ) 

where 

5 _ (A + 2 G) 

" 1 - A Hi ’ 

and A and /xx are the area and area moment of inertia about the x axis, respectively. 
Similarly, the local shear resultant can be expressed as 

Q(z,/) = | TyzdA = koG A{v lZ +9) , (5.a) 


I 


= I y (U 0 , x + V 0<y ) dA 


(4 .b,c) 


where, M,z * Q, 

1(0 = (1 -gh 2 ) 


(5 .b) 
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and 

H 2 « ^(V 0 +W 0 ,y)dA. (5.c) 

Substituting Eqns. (4 .a) and (5.a) into Eqn. (1 ) results in the final form of the displacement 
relations defined in terms of the kinematic variables (v and 6) only; 

u(x,y,z,f) = 6{z,t), z Yx{x,y), 

v(x,y,z,t) = v{z,f) + eiz,t), z \f/y{x,y), (6 .at-c) 

w(x,y,z,f) » y 9(z,1) + [v{z,t), z +6{z,())yr z (x,y) 

where 

Yx(x,y) » El U 0 (x,y) , 

Yyix,tf = El V 0 {x,y ) , (6-d-O 

Yz(x<y) - k 0 G A Wo(x,y) . 

The six strain components can be written by substituting Eqns. (6.a-c) into (2.a-f); 

£ xx = 9< Z Yx,x . 

£ yy - Q’Z Yy.y . 

£ zz = yO,z + {v, z + e), z y z , 

(7.a-f) 

7yz - { v <z + 0)(1 + Yz.y) + &zz Vy » 

7xz = (v.z + ^XV'z.x) + #>zz Yx . 

7xy = Qfd^Yx.y + V^y.x) • 

The equations of motion are derived via Hamilton's principle 

f'2 

<5/7 » (<5(/-<5T-<5W e )cfr = 0 (8) 

■"l 

where SU, ST, and SWq are the variations of the strain energy, the kinetic energy, and 
the work of external forces, respectively. The variation of the strain energy, which can 
be written in two parts, is given as; 
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su 


r \^ T[c] 


dA dz ■ 5Ub + SUs 


(9. a) 


where the first portion represents the strain energy variation associated with bending; 


SU B = 




1 se, 7 Ffo,, o, 2 


1 Sv, zz 1 1 01 2 D 22 . 

v.zz 

Jo 

^ J 


dz, 


(9 -b) 


the second portion is associated with shear; 


5U S = 



5v, z + 59 

59, zz 


7T 


^11 2 1/ 

A\z A22 Jl 


v, z + 9 
9,zz 



(9.c) 


and the coefficients Ajj and Djj (i,j=1,2) are defined in the Appendix (A. 1-6). The 
constants A-\ 2 and A 22 represent additional shear deformation as a result of in-plane 
cross-section deformation. The variation of the kinetic energy is given as 


5T 



Su-u + 5v‘V+5w-w\dAdz 


•vvjcf/ 1 


( 10 ) 


where p is the mass density and ( ) represents a derivative with respect to time; t. 
Substituting Eqns. (6.a-c) into (10), taking the appropriate derivatives with respect to 
time; t, and then integrating over the cross-section; A, the resulting kinetic energy 
variation is equal to; 


5T 


r - 

5v ' 

T 

i 

8\t, z 

1 


59 


J 

j 89, z 

- 


m 
0 
0 
J 1 


0 
Jz 
J 2 
0 


0 


Jz 

b 

0 


Ji 

0 

0 

j p 


V 


v'.z 

9 

e.z 


)dz 


1 


( 11 ) 
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and the nonzero matrix coefficients are defined in the Appendix (A.7-12). The constants 
Ji , and Jp represent the additional kinetic energy associated with in-plane cross- 
section deformation. The variation of the work of external forces is given as; 


( a*. 

Jo 


5W e » I p(z,tjSv a dz, 


(12) 


where va represents the displacement on the surface of the beam at the point of the 
applied load ( v(x=0 ,y=y,z,t)). 

The two differential equations of motion and associated boundary conditions are 
obtained by substituting Eqns. (9.a-c), (10), and (12) into (8) and integrating by parts; 


(Ai •)(/,/+ 0) + Ai20,zz},7 ■ [D22^>zz + 0i20>z}>zz- mV + J-\9,z 


JzV.z + JzQ 


>z - P > 
(13.a) 


(Ai i (v,z + 0) + Ai20.zz} "(011 0.z + D-\2 v <zz)>z + (A220.zz +■ A 12 ( v <z + 0)}»zz 


|ji V + Jp0,zj,z ’ P<2 


(13 .b) 


= -/p0 - J 2 V.Z + Ji V + Jp9,z\z ■ P<z¥y 


where VV is the magnitude of the warping displacement at the point of the applied load 
given as 

¥y = Vy(x=0,y=y) . (13.c) 

There are four sets of boundary conditions that must be specified at the beam ends 
(x=0,L). First, define either the transverse displacement (v) or the effective shear force 

Q » (Ai \(v, z + 0) + A|20.zz} ' (D 22 V'.zz + 0i20«z}«z 


+ Jz^.Z + J2S • 


(14.a) 


Second, specify either the rotation of the cross section (0) or the effective cross-section 
moment; 


- {^220.zz + A-\ 2 {v, z + 0)}.z" PVy • 


(14 .b) 


Third, define either the derivative of the transverse displacement (v,*) or the generalized 
moment; 

M — Doqv .77 +• Pi (14.c) 

Finally, specify either the derivative of the rotation {$,£) or; 

R s A22&izz + Ai2(v,z + 6) . (14. d) 


3. SIMPLIFIED GENERAL HIGHER-ORDER THEORY 


A simplified version of the general higher-order theory can be developed by 
assuming that the cross-section is rigid within it's own plane (i.e. f the third displacement 
group of Eq. (1) is zero ( U 0 , V'o-O)) and that the stresses within the arbitrary cross-section 
are negligible (axx.Cyy.Txy-O). This model represents a logical extension to existing 
higher-order theories [11-14], which are limited to thin rectangular cross-sections. The 
displacement field will have the form; 

u{x,y,z,t ) » 0 , 

v(x,y,z,t) = v{z,t ) , (15.a-c) 

vv(x,y,z,f) = y 9{z,t) + (v(z,/), z +0(z,f))v'z(*,y) 


where 

and 


V'ztoy) = k 0 G A W 0 {x,y) , 


ko — 



Wo,ydA 



(15 .d) 
(15.e) 


The resulting three nonzero strain components have the form; 


(16.a-c) 


ezz - yO,z + (v,z + 0\z Yz , 

Yyz - {v,z + 0)(1 + Yz.y ) . 
yxz - {v,z + &fVz,x) • 

The one-dimensional constitutive relations are developed by assuming that the stresses 
within the cross-section are zero (c^x,Oyy,r xy - 0 ) in Eq. (3.a), then a static condensation 
approach is used to solve for the remaining three stresses in terms of the corresponding 
three strains; 




E 0 0 

&zz \ 

*** 

l 1 

| 3 

0 G 0 
.0 0 G. 

Yyz 

Yxz 1 


(16 .d) 


where E and G are Young's modulus and the shear modulus, respectively. 

The equations of motion are again derived via Hamilton’s principle (Eq. (8)), where 
the simplified form of the strain energy variation is expressed as; 


5U = 


I 59,z j 7 " 

Du D12 ( 

M 

1 5 v,zz 1 
yo 

. D12 D22J 

v.zz j 


+ (8v>z + 8Q)T Ay i(v,2 + 9)dz , (17. a) 


the simplified form of the kinetic energy is given as; 



f L 

1 

s* 

T 

If 

*o 

( 

S*.z 

\ 
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j 

0 

1 



m 0 0 

0 Jz Jz 

0 J 2 Ip 


v 

fr.z )<*z 

9 


(17 .b) 


the simplified form of the work of external forces is defined as; 


•r 


SW Q = I p(z,f)8vdz , 

Jo 


(17 C) 


and the associated section constants are given in the Appendix (A. 13-20). 
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The resulting differential equations of motion have the form; 

{Aii(v,* + 0 )}, 2 * [D22 v >zz + 0l20.z}.zz = nN~ |«^zV»z + ("1 9. a) 

A i^z + 9) * { 0-| 1 9,z + D-\2V>zz)’Z m m lffl~'hy*z i (18.£>) 

where the three geometric and natural boundary conditions that must be specified at the 
beam ends (x^O ,L) include; 

fiflQnrct ns Naluiai 

v Q = /Aii(v, 7+ 0)-|D22' i ',zr + Oi28,zj,z+ JzV,^+ ^20 

v,z M =* Dzz v >zz + ^I20»z (18.c) 

Q M = Di 1 9,z + D-\2 v >zz 


4. MODEL VERIFICATION 


The differential equations of motion for the general higher-order theory (Eqns. 13.a,b) 
and the simplified general higher-order theory (Eqns. 18.a,b) reduce identically to the 
linear form of the higher-order theory for beams, developed by Heyliger and Reddy [14], 
if one assumes a rigid in-plane cross-section {yx = yy = 0), negligible in-plane stresses 
(oxx,Oyy.'rxy*0), and a cross-section that is an extremely thin (plane stress) rectangle so 
that the out-of-plane shear-dependent warping can be described from Saint-Venant's 
static flexure warping function [23] as; 


Yz = 



09. a) 


where (2 b) is the overall height of the cross-section (Fig. 2) and the remaining eight 
section constants are determined by making use of Eqns. (A13-20); 


m = pA, 


Au - ^GA , 


On 


= -£S-e/ xx , 

105 


012 
0 22 


^El xx , 

105 

^-El xx , 
21 


Ip « 
Jz 
Jz 


plxx , 
1 05^ 

fo! p '"' 

1 -- plxx • 

2T 


(19.6-/) 


I 


Furthermore, the current equations will reduce to the well-known Timoshenko equations 
of motion [1,2] by neglecting both in-plane deformation {y* = Vy * 0) and out-of-plane 


warping {y z = 0) effects; 


■ 

[A-\-\(v,z+ 0)],/ ■ mV- p , 

(20. a) 




■ 

(Ai i(v, z + 0)] - (Du 9,z)tz - > 

(20.6) 



where the four section constant are equal to: 


1 

An = kGA , m = pA, 

0i i * Elxx i Ip - plxx > 

(20. c-/) 

1 

with k being a shear coefficient that is dependent upon the geometry and 

material 

1 


definition of the cross-section [1-7]. 


5. FREE VIBRATION OF A SIMPLY-SUPPORTED BEAM 

The natural frequencies of a simply-supported beam can be obtained by applying the 
Ritz method where the following displacement functions satisfy both the geometric and 
natural boundary conditions; 

v{z,fl = Qi sinJn^Jsin(o?f), 9{z,t) = Qz cosjn^j sin(cof), (21. a, 6) 

where n is equal to the mode number, co is the corresponding natural frequency, and Qi 
and 02 are unkown generalized coordinates. Alternatively, this solution represents the 
motion of a standing wave in an infinitely long beam where n=1 and L is taken as the 
wavelength. Clearly, as L becomes small (or n becomes very large), the effects of local 
in-plane deformation and out-of-plane warping (shear-deformation) will become 
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significant. Substituting the assumed displacement field (21.a,b) into the strain energy 
(9.a-c) and kinetic energy variations (11), integrate over the length; L, and substituting 
the results into Hamilton's principle (8), will lead to the following 2 by 2 matrix form of the 
algebraic eigenvalue equation; 


([K]-«2[M]}(0)-(0| 


(22.a) 


where, the coefficients of [K\ are given as 


with 


*ii = £ 




*12 = *21 = ^ 


*22 * ^ j ^ 11 + ^ Dl1 ’ 2 > 4 1 2 ^) + ^ 22 (^ f ] 


Wi- 


ll 

2 


m+ J; 


w 

(Jz ~ j,H 


(j 2 - 

iff 


In + Jr 


(22 .b) 


(22. c) 


and 


(Q| r -(Q 1t Oz ) . 


(22 .a) 


The coefficients of the stiffness and mass matrices for the simplified model can be 
obtained by setting (A-\2, A 22 . >^ 1 . Jp) equal to zero in (Eqns. 21.b,c) and calculating the 
remaining eight section constants using (Eqns. A. 13-20). 


6. NUMERICAL EXAMPLES 

Two sets of numerical examples are presented. First, the general higher-order 
theory (GHOT) and the simplified higher-order theory (SGHOT) are validated with 
previously published solutions over a broad range of beam (or wave) lengths. Second, 
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the vibration behavior of a simply-supported beam is studied for a wide range of 
rectangular and elliptical cross-section aspect ratios (b/a). See Fig. 2 .a,b. 
Comparisons with existing models are made so that the features of the current 
development can be assessed. The required section constants for the rectangular and 
elliptical cross-sections are presented in tables 1 and 2, respectively, where the in- 
plane deformation functions (y/*, y/y) and out-of-plane warping function (yr z ) were 
developed from Sokolnikoffs full three-dimensional elasticity solutions [23], assuming 
that the Poisson's ratio (v) is equal 0.333 and b » 1. For example, the force-dependent 
warping functions for a rectangular cross-section have the form 


U 0 {x,/) =--=£- xy, 
E Ixx 


(23 .a) 


Vo{x ' y)s ~TE]7x {y2 ' x2) ' 


(23 .b) 


Wc 


(x.y) * wh l^y*-fyx2 til 

c Ixx 6 2 K 3 n 3 

n- 1 " 


Z, (a sinh&L) cos(Pz*-)-nny) 


_a 

cosh 


, (23. c) 


and the constant ( kp ) is calculated using Eq. (5.b) as; 


kc — 


1 + 


1 


2 (1 +v)' 


1+v(1- f ) 


Jafi-iajkfv tnl 1 


Vi 


INI " 2 


costing 



(23. d) 


This coefficient (ko) will approach 2/3 for very thin rectangular cross sections (a/b» 0). 
The desired form of the warping functions [yx,¥yt¥z) can be determined by substituting 
Eqns. (23 .a-d) into (6 .d-f). Thus, 


Yx = ~vxy , 


(24 .a) 


VO^-jr-O' 2 '* 2 ) • 


(24.b) 


and ¥z is found using Eq. (6./), where Wo and kp are given in Eqns. (23 .c,d). 
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In the following figures, the two calculated frequencies (o>i, ©2 ) are normalized 
with respect to oq (= m(^Bp)/L) and plotted versus (2 nb/L). For extremely long slender 
beams (2nb/L-0), the lower frequency, co -\ , is just the Bernoulli-Euler prediction with 
(Qi*0, 02=0) , while the upper frequency, o>2. will equal the "thickness-shear" 
frequency, cats (=VAii//p), where the beam experiences pure shearing through the 
cross-section (02*0) with no beam bending (Gi=0). As (2 nb/L) increases, the modes 
become coupled, where the lower natural frequency, ©i , is predominantly beam 
bending with some shear deformation, whereas the upper frequency, g> 2, is 
predominantly beam shear with some beam bending. A second (lower) x-axis, is 
included in the figures, that represents the corresponding mode number assuming a 
characteristic long slender beam (2b/L=0.1). 

6.1 VALIDATION STUDIES 

The current theories; (GHOT) and (SGHOT), were validated by comparing the 
calculated natural frequencies of a simply-supported beam having a thin rectangular 
(b/a=1000) cross-section with numerous existing solutions over a broad range of beam 
lengths (Fig. 3) and by comparing the calculated "thickness-shear" frequencies (Table 
3) with previously published results [24]. The lower frequency co i , has four distinct lines 
which represent; (1) the Bernoulli-Euler solution, (2) the current (GHOT) solution, (3) a 
plane-stress elasticity solution [7], and (4) the Timoshenko solution with fc=2/3 [1,2]. The 
line which contains the plane stress elasticity solution also represents the current 
(SGHOT) model, a higher-order theory [14], and the Timoshenko solution with k =. 8551 
[6]. All of the shear-deformable models are in complete agreement with the Bernoulli- 
Euler solution when the beam is long and slender, but as the beam (or wave) length gets 
shorter, the Timoshenko beam theory with k= 2/3 is more flexible than the remaining beam 
theories, and for extremely short beam lengths (2b>L) the current (GHOT) model predicts 
a natural frequency that is slightly higher than existing one-dimensional beam theories. 
This difference is clearly a result of including in-plane deformations in both the strain and 
kinetic energies. 

The shear-dominant frequency eo2, has only two distinct lines which represent; (1) 
the Timoshenko solution with k=2/3 [1,2], and (2) the Timoshenko solution with fc=. 8551 
[6]. The second line also represents a higher-order theory [14], and the current (GHOT) 
and (SGHOT) models. The Timoshenko solution with k=2/3 is measurably more flexible 
when the beam is long and slender, but converges to the other beam solutions when the 
beam (or wave) length becomes very short. 
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A comparison of the predicted "thickness-shear" mode frequency (Table 3) with 
an elasticity solution [24] shows that the current (SGHOT) model and the model of [1 1] is 
in near exact agreement, while the Timoshenko-based predictions are somewhat 
deficient. The current (GHOT) model differs slightly as a result of including the kinetic 
energy associated with in-plane deformations, which was ignored in [24]. 


6.2 VARYING CROSS-SECTION ASPECT 

In the second study, the vibration behavior of a simply-supported beam was 
analyzed for a wide range of rectangular and elliptical cross-section aspect ratios ( b/a ). 
In Fig. 4, the two calculated natural frequencies using the current (GHOT) model are 
presented as a function of rectangular cross-section aspect ratio (b/a). It can be seen 
that (b/a) has very little effect on o>i , even for extremely short beam (or wave) lengths, 
but it can have a measurable effect on o>2 for moderate and shorter beam lengths. In 
Table 4, a comparison of the lower and higher frequencies is given as a function of b/a 
for an extremely short beam length (2 b=L) and one can see that b/a has a minimal effect 
on the bending frequency, whereas the effect on the shear-dominant frequency can be 
significant. 

As a further validation of the current model, Aalami [15] predicted that ar\ =0.4724<ao for 
an extremely short beam (2 b=L) having a square cross-section (b/a= 1) using a 
Rayleigh-Ritz energy approach on the full three-dimensional problem, which is within 
0.4% of the current (GHOT) solution. The "thickness-shear" mode frequency is also 
presented in Table 4, where it is seen that the frequency increases with a reduction in 
the aspect ratio (wider cross-section). The results from the Bernoulli-Euler solution, the 
Timoshenko solutions [1,2,6], and the higher-order models [11-14] were not included in 
either Fig. 4 or Table 4 because these solutions do not depend upon aspect-ratio for 
rectangular cross-sections and are equal to the results given in Fig. 3 and Table 3, 
respectively. 

In Figs. 5-7, the two calculated natural frequencies are presented as a function of 
elliptical cross-section aspect ratio (b/a) for the Timoshenko model and the current 
(GHOT) and (SGHOT) models. The shear correction factor (k) used in the Timoshenko 
model was taken from [6], where the magnitude varies with aspect ratio. The results from 
the Timoshenko solution (Rg. 5) show that ar\ and a>2 become slightly more flexible with 
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decreasing aspect ratio (b/a), where a large unrealistic reduction occurs when 
b/a< 0.50. This reduction is apparent in short wave-lengths for the bending-dominant 
frequency (<wi) and in the long wave-lengths for the shear-dominant frequency (a>2). 
Alternatively, the frequency results using the current (GHOT) model (Fig. 6) show very 
little change in ©i with b/a except for extremely short beams (2 b>L), whereas the 
Change in co2 can be significant over the broad range of wave-lengths. It is interesting to 
note that aspect ratio has a negligible effect on both the calculated frequencies for a 
wide range of cross-section shapes (1<b/a<1000), but as the cross-section becomes 
wide (b/a<1.0) the variation in the calculated frequencies, especially co 2, can be 
important. Finally, the frequencies as a function of beam-length and aspect ratio b/a 
using the current (SGHOT) model are presented in Fig. 7, where the bending dominant 
frequency (o>i) will undergo a slight increase for extremely short beams ( 2b/L>2 ), while 
the shear dominant frequency (co2) experiences a small decrease for long slender 
beams. 

In Tables 5.a,b, g>i and <02 for an extremely short (2b=L) simply-supported beam 
with an elliptical cross-section are compared using the three different solution 
approaches for a variety of aspect ratios. It is clear that decreasing b/a (wide section) 
will decrease both 01 and C 02 based upon either the Timoshenko results or the (SGHOT) 
model, but the results using the (GHOT) model show that the co 1 may actually increase 
for wide cross-sections (low b/a), which is in agreement with the three-dimensional 
results using Aalami's model [15]. From Fig. 6, it is apparent that as the beam wave- 
length is much smaller than the cross-section dimensions (2b), then (< 01 ) will decrease. 
Finally, in Table 5.c the "thickness-shear" mode frequency is presented using the three 
different methods and compared to an elasticity solution of Jeffreys [25]. Both the 
(GHOT) and the (SGHOT) models are in near perfect agreement with the exact results for 
a wide variety of cross-section aspect ratios, whereas the Timoshenko based solutions 
[6] slightly over predict the results. 


7. CONCLUSIONS 

A new one-dimensional theory has been developed to study the vibrational 
behavior of prismatic beam-type structures that utilizes both out-of-plane shear- 
dependent warping and in-plane (anticlastic) deformations. The two equations of motion 
were derived using Hamilton's principle, where the full three-dimensional constitutive 
equations are used and it is assumed that the cross-section deformation functions are 
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Known. A simplified form of the current general higher-order theory was also presented 
which includes out-of-plane shear-deformation but neglects In-plane deformations by 
assuming that the stresses within the cross-section were negligible. It was shown that 
the current model reduces identically to existing higher-order models of beams having a 
thin rectangular cross-section, when in-plane deformations are not included. 
Furthermore, the current model reduces to the Timoshenko equations when out-of-plane 
shear-dependent warping is neglected. Results from a numerical validation study of a 
simply-supported beam with a thin- rectangular (b/a=1000) proved that the current model 
is in near exact agreement with existing approaches over a broad range of beam (or 
wave) lengths for the bending-dominant, shear-dominant, and thickness-shear 
frequencies. It was shown that including the in-plane cross-section deformations and the 
three-dimensional constitutive model will slightly increase the bending-dominant 
frequency and greatly reduce the shear-dominant frequency for extremely short wave- 
lengths. Results from a second numerical study showed that the cross-section aspect 
ratio has only a minimal effect on the bending-dominant frequency even for short wave- 
lengths, but the effect on the shear-dominant frequency can be significant for long 
slender beams with either rectangular and elliptical cross-sections. Moreover, the 
"thickness-shear" mode frequency was found to increase slightly with decreasing aspect 
ratio for the rectangular cross-section, whereas the frequency associated with this mode 
was found to decrease for the elliptical cross-section and be in near exact agreement 
with existing published solutions. 
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APPENDIX 


\ 


The bending and shear related cross-section constants for the current general higher- 


— 

order theory (GHOT) are defined as 



— 

On = | (A+2 G)|(y v^ x + v^ t/ J+2 aJ(/ + yr 2 | y Kx + v/y, y )+ V'** V%,)+g{ Yk y +Vy, x j 
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The mass related cross-section constants for the general higher-order theory (GHOT) 


are defined as 
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and 


Jp * | (^Yx + vf) CM 

(A 10) 

Jz = f PY'f 

(All) 
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^2 = | p(y+ Yz)Yz dA . 

(A 12) 


The following eight nonzero section constants are used in the simplified general higher- 
order theory (SGHOT) 
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Table 1 Nondimensionalized section constants for rectangular cross-sections (b/a) 
with v* 0.333. 


b/a 

0.50 

1 

2 

10 

100 

1000 

A<\\IGA 

.55373 

.60373 

.58046 

.57019 

.56977 

.56976 

A]2lGAk£ 

.10788 

.0075697 

-.018517 

-.026057 

-.026356 

-.026359 

AvilGAbfi 

.11890 

.017249 

.007431 1 

.0056066 

.0055451 

.0055445 

Oll/E/xx 

.67878 

.73949 

.71522 

.70501 

.70460 

.70459 


-.12390 

-.10703 

-.11659 

-.12024 

-.12039 

-.12039 

D22^hix 

.07343 

.046448 

.051600 

.054504 

.054630 

.054631 

m/(pA) 

1 . 

1 . 

1 . 

1 . 

1 . 

1 . 

/p/(pfxx) 

.65436 

.72405 

.69806 

.68689 

.68643 

.68643 

^1 /(p5cx) 

.49950 

0 

-.12488 

-.16484 

-.16648 

-.16650 

JlUph cx) 

-.14831 

-.12248 

-.13375 

-.13837 

-.13855 

-.13855 

Jzf(phtx) 

.049015 

.031004 

.034443 

.036381 

.036465 

.036466 

Jp/iph cx^) 

.35669 

.051748 

.022293 

.016820 

.016635 

.016633 
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Table 2 Nondimensionalized section constants for elliptical cross-sections (b/a) with v 
m 0.333 for (a.) the general higher-order theory, and (b.) the simplified general 
higher-order theory (where A-\ 2 » A22 * ^1 = Jp * 0.0). 


b/a 

0.50 

1 

2 

10 

100 

1000 

A^IGA 

.41945 

.57846 

.64773 

.67502 

.67622 

.67654 

A\ 2/GAb^ 

.075454 

.015139 

-.013469 

-.024242 

-.024707 

-.024712 

A22)^‘Abfi t 

.068151 

.0092408 

.0042599 

.0034887 

.0034655 

.0034653 

Dll/E/xx 

.58259 

.69203 

.74236 

.76257 

.76345 

.76346 

Di 2/Ekx 

-.12992 

-.12612 

-.10940 

-.10146 

-.10109 

-.10109 

D22/EIxx 

.15757 

.055736 

.038845 

.034518 

.034354 

.034352 

m/{pA) 

1 . 

1 . 

1 . 

1 . 

1 . 

1 . 

tyiplxx) 

.53020 

.67350 

.72944 

.75109 

.75203 

.75204 

W(pkx) 

.49950 

0.0 

-.12488 

-.16483 

-.16648 

-.16650 

J2/{plxx) 

-.18231 

-.14465 

-.12231 

-.11293 

-.11251 

-.11251 

JzKpIxx) 

.10518 

.037204 

.025929 

.023041 

.022931 

.022930 

J^plxxt^) 

.27260 

.036963 

.017038 

.013955 

.013862 

.013861 


b/a 

0.50 

1 

2 

10 

100 

1000 

A-\ i/G/4 

.42130 

.51852 

.58848 

.62334 

.62498 


011/E/xx 

.53472 

.61111 

.66975 

.69994 

.70137 

.70138 

Di2/E/xx 

-.18750 

-.16667 

-.14506 

-.13257 

-.13195 

-.13195 

022 lEIxx 

.090278 

.055555 

.040123 

.034931 

.034724 

.034722 

m /( pA ) 

1 . 

1 . 

1 . 

1 . 

1 . 

1 . 

tyiplxx) 

.53472 

.61111 

.66975 

.69994 

.70137 

.70138 

J2t{pkx) 

-.18750 

-.16667 

-.14506 

-.13257 

-.13195 

-.13195 

Jz/(plxx) 

.090278 

.055555 

.040123 

.034931 

.034724 

.034722 
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Table 3 Natural frequencies of a short (2 b/L * 1) simply-supported beam with a thin 
rectangular cross-section (b/a = 1000) including thickness-shear mode 
frequency (cots)- 



elasticity 

Bernoulli- 

Euler 

Timoshenko [1 ,2] 
k = 2/3 k « .8511 

Refined 

111]- 

KsUEjI 

■gSSSSI 

oy\/coo 

.4666 [7] 

.9069 

.4269 

.4623 

.4628 

.4763 

.4628 

g/2./(Dq 

- 

- 

.1712 

1.2222 

1.2138 

1.2010 

1.2138 


.9069 [20] 

- 

.8166 

.9226 

.9076 

.9111 

.9076 






Table 4 Natural frequencies of a simply-supported beam (2 b/L * 1) for various 
rectangular cross-section aspect ratios. 



b/a = 1/2 b/a= 1 tfa = 2 #a=1000 

© 1 / 0)0 
qjz/qjq 
0 ) ts / t»S 

.4786 .4744 .4756 .4763 

.9762 1.1486 1.1890 1.201 

.9200 .9132 .9120 .9111 




Table 5 Natural frequencies of a simply-supported beam (2b/ L » 1) for various elliptical 
cross-section aspect ratios for (a.) the bending dominant frequency (o^/oio), 
(b.) the shear-dominant frequency (coz/(Oo), and (c.) the thickness-shear 
frequency (^s/cos). 



Bernoulli- 

Euler 

Timoshenko 

[6] 

mmSSM 

WmmSSBm 

bla » 0.50 

.7854 

.4377 

.4636 

.4394 

iVa» 1 

.7854 

.4462 

.4563 

.4436 

b/a = 2 

.7854 

.4489 

.4564 

.4474 

b/a= 1000 

.7854 

.4499 

.4566 

.4495 



Timoshenko 

[6] 

Current 

(GHOT) 

Current 

(SGHOT) 

b/a = 0.50 

1.274 

1.014 

■■BK • ^maam 

b/a = 1 

1.294 

1.209 


b/a = 2 

1.301 

1.264 

1.290 

b/a = 1000 

1.304 

1.282 

1.294 



elasticity 

[3,21] 

Timoshenko 

n 

■m 

WmssSSSM 

b/a = 0.50 

1.257 

1.289 

1.257 

1.255 

b/a = 1 

1.303 

1.333 

1.305 

1.303 

b/a = 2 

1.328 

1.349 

1.331 

1.335 

b/a = 1000 

1.339 

1.354 

1.342 

1.336 










(x = 0 ,y=y,z) 



Fig. 1 A prismatic beam. 
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mode number (2b/L = 0.1) 


Fig. 3 Natural frequencies of a simply-supported beam with a thin rectangular cross 

section (b/a=1000) ( — - — [6,7,11] and SGHOT, GHOT, [1,2], — 

Bernoulli-Euler). 
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Fig. 4 Natural frequencies of a simply-supported beam with a rectangular cross- 

section using the general higher-order theory (GHOT), ( #a=10-1000, 

b/a= 2, tfa=1 , — b/a= 0.5, Bemoulli-Euler). 


532. 
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« i l i i » i ; • j 

1 2 10 

mode number (2b/L = 0.1) 


40 


Fig. 5 Natural frequencies of a simply-supported beam with an elliptical cross-section 

using the Cowper's version of the Timoshenko theory [6] ( b/a=1000 

(k=0.91718), b/a=.S (k= 0.8296), b/a=0A {k= 0.3040), 

Bernoulli-Euler). 
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Fig. 6 Natural frequencies of a simply-supported beam with an elliptical cross-section 

using the general higher-order theory (GHOT), ( iVasl-IGOO, 

b/a=.5, — Bernoulli-Euler). ^ : • 




CO/CD 



mode number (2b/L = 0.1) 


gg Fig. 7 Natural frequencies of a simply-supported beam with an elliptical cross-section 

using the simplified general higher-order theory (SGHOT), ( b/a=1-1000, 

g b/a=.5, Bemoulli-Euler). 
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Chapter 6 

FORMULATION OF A NONLINEAR THEORY FOR SPINNING ANISOTROPIC BEAMS 

BY 

C.A. Ie and J.B. Kosmatka 

Department of Applied Mechanics and Engineering Sciences 
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ABS TR AC T 


A geometrically nonlinear theory is developed for spinning anisotropic beams having arbitrary 
cross-sections. An assumed displacement field is developed using the standard three-dimensional 
kinematic relations to describe the global beam behavior supplemented with an additional field that 
represents the local deformation wi thin the cross section and warping out of the cross section plane. It 
is assumed that the magnitude of this additional field is directly proportional to the local stress 
resultants.Using a developed ordering scheme, the nonlinear strains are calculated to the third order. 
Through Hamilton's principle, the six governing equations are obtained.The finite element model is 
developed using the weak form variational formulation. Numerical results for a static case show that 
the model agrees with the elasticity solution up to the stress level. Results on the free vibration cases 
show that the behavior of am so tropic beams are indeed more complex compared to the isotropic 
counterpart, i.e. complete coupling (bending torsion shear and extension) exists for these type of 
beams. 


INTRODUCTION 


Consider a prismatic beam of length (L) composed of an arbitrary anisotropic material having a 
general cross-section (see Fig. 1). For solely the purpose of simplicity of formulation, but without 
loss of generality, let us assume that the material is homogeneous through out the entire body 
Cartesian coordinate systems (X, Y, Z) and (x, y, z) are defined on the beam where y and z are 
coincident with the principle axes of the root cross-section and x is coincident with the line of 
centroids. The beam is located (//*, Hy, H z ) with respect to the coordinate system (X, Y, Z) and 

spins about the coordinate Z with a constant angular velocity Q. The beam is subjected to distributed 
loads px, py, and p z , which are assumed to act on the coordinate line; x . 

Existing beam theories have been developed by neglecting the in-plane stresses, i.e., <Syy, Ozz 

c yz . (Love, 1927). Classical beam theory further assumes that plane sections perpendicular to 
the undeformed x axis remains planes (first assumption) and perpendicular (second assumption) to the 
the deformed x axis. As a consequent, the shear strains will automatically vanish. To take into 
account the shear strain energy, Timoshenko (1921, 1922) developed a theory that abundant the 



second assumption. Due to the fact that these theories were developed for isotropic materials, they 
assume uncoupled global material behavior. 

Beams theories that were developed by abandoning both constraints exist Levinson (1981), 
Heyliger and Reddy (1988) and Kant and Gupta (1988) developed theories that more accurately 
represent the kinematical relations by introducing shear related warping function which is proportional 
to the intensity of the shear strains at the centroid These theories were developed exclusively for 
isotropic material with a thin rectangular cross section. Vlasov (1961) developed a theory that also 
include the out-of-plane warping for thin-walled isotropic beams having simple cross-sections. 
Bauchau (1985) extended this approach for thin- wall composite beams where eigen- warping functions 
were used to model the out-of-plane shear-dependent warping. Recently, Ie and Kosmatka (1992.a) 
developed a beam theory for isotropic general cross-sections using first-order warping functions that 
can predict not only the global behaviors but also the local behaviors (strains and stresses). 

This paper extends the development by Ie and Kosmatka (1992.a) to a more complex material 
definition, i.e. anisotropic materials. For the sake of simplicity and clarity, but without loss of 
generality, this theory is developed for the case of beams having homogeneous material through out 
the body. A set of numerical examples is presented for beams having elliptical cross sections. Exact 
warping functions for this cross section can be found in a paper by Ie and Kosmatka (1992.b). For 
beams having general cross-sections and material definitions, the warping functions can be found 
using a model developed by Kosmatka and Dong (1991). 



Fig. 1 A spinning anisotropic beam. 


THEORETICAL DEVELOPMENT 


Assumed Displacemen t Field. Mix«d Formulation 


Through out the body, the displacement distributions in the x -, y- and z -directions are defined as 

6 . 2 - 


u(x,y,z,t)=u(x,t)+z Q y (x,t)-y 0 Z (JC,I)+W 1 (xo>,z t i) , 
v (x,y,z,t) =v(x,t)-z Q x (x,f) + W 2 (x,y,z,t) , 
w ( x,y,z,t ) =w(x,t)+y 0 x (jr,r) +W 3 (x,y,z,r) , 


(l.a-c) 


where u , v, and w represent the displacements in the directions of the coordinate axes x-, y-, and z- 

respectively and Q x , 0 y and 0 2 represent the rotations about the coordinate axes x-, y-, and z- 
respectively. Wj,W2 and Wj, represent the total local warping functions (displacements) of the 
cross sections in the directions of die coordinate axes x-, y-, and z- respectively . These functions are 
assumed to have the following form 


where 


6 

Wj (x,y,z,r) = £ Fj(x,r) w u (y,z) , 
1*1 

6 

w 2 ( x,y,z,t ) = £ Fi (x,t) w 2 i (y,z) , 

i =1 
6 


W 3 (x,y,z,t) = £ Fi (x,t) w 3 , (y,z) , 



i =1 

Fi(x,t) = I Cxx dA , 

F 2 (x,t) = I o xy dA , 

JA 

f ' 

JA 

F 3 (x,r)= I o zx dA , 

F*(x,t) =1 {y a zx - z o xy ) dA , 

JA 

r 

JA 

* 

Fs(x,f) =1 zOxxdA , 

F 6 (x,r) = - I yo xx dA , 

JA 

JA 


(2.o-c) 


(3 .a-f) 


in which A is the area of the cross section, Oxx, Qyy and a^z 3Te the normal stress components and 

Oyz, Ozx 311(1 a xy are the shear stress components. Fj , F2 and Fj , are the local stress resultants 
(forces) in the x-, y-, and z- directions respectively and F4 , F5 and F5 , are the local stress resultants 
(moments) about the coordinate axes x, y, and z respectively. Further, wj( , W2/ and w$j (named 
as beam warping functions) are the warpings of the cross section with a thickness dx (zero thickness) 
due to individual stress resultant F\ provided that the area dA is fixed against translational 
displacements and rotations. Solving the boundary value problem given in Fig. 2 , in which the only 
nonzero stress resultant is Fj, the warping functions w/j , W2J and wjj will be obtained. The 
others fifteen warping functions can be obtained in a similar way. 



Fig. 2 A procedure to obtain warping functions. 
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The general constitutive relations are given as 


{a}-[C](e) f 


(4) 


where the stress and strain arrays are defined as 

{o) r = { Oxx , Gyy, Gzz, Gyz, Gzx, <Jxy } , 
{e } T — { £xx, tyy,tzz, \yz, Y zx, yxy } , 


(5 Ab) 


and [C] is the material stiffness matrix which may have twenty one distinct stiffness constants. 


Assumed Displacement Field Purely Kinematic Formulation 

It is necessary to transform the assumed displacement field in eqns (l.a-c) to a purely kinematic 
field counterpart before further formulations can be furnished. For this purpose, a linear strain- 
displacement relationship is being adopted, Le. 


du 

dx 

dw . dv 






_ dv 

£yy ~dy ^ 
du , dw 

Y ~ = al + e^ 


_ dw 
dv du 


Yyz dy dz ' 

On the basis of beam elemental equilibrium, the following assumptions are being made 


(6.a-/) 


Fj — Px — 0 1 F 2 — Py ~ ® > 

F 3 = -p z = 0 , F 4 = 0 . 

Also note that 


(7 Mril) 


F 5 = F 3 , F 6 = -F 2 . (J.eJ) 

Using the above assumptions, the local stress resultants can be expressed in terms of u, v, w ,Q X , Qy 

and Q z (and their derivatives) by first calculating the strains using eqn (6 m-J) and eqn (l.o-c), then 
using this strains to get the stresses through eqn (4), and finally substituting the calculated stresses into 
eqns (3 .a-J) to get the following result 


where{F } is an array 


{F)= [H}{e), 


( 8 ) 


(F)7'={F 1 ,F 2 ,F 3 ,F4,F5,F 6 ) , (9) 

[H\ is a six by six matrix which depends on the geometry of the cross section, material properties 
and warping functions, and (e } is an array 

{^ = { ^l» ^2» ^4* } (10) 

= («’, V-0 Z , w' + 0y, 0 X , 0y, 0 Z ) . 

Finally, the desired form of the assumed displacement field is obtained by substituting eqn (8) into eqn 
(1) through eqn (2), resulting in the following formulation 


6. 4 



u ( x,y,z,t ) = u(x,t)+z Qy(x,t)-y 0 z (x,r)+ipi(x,y,z,r) , 
v( ,y,z,t) = v(x,f)-z 0*(x,r) +t|> 2 (x, y,z,f) , 
w(,y,z,t) = Mx,t)+y 0 x(x,r)+H> 3 (j:,y,z,f) , 


(1 1 .a-c) 


where 


6 

V. (x,y,z,t) = £ e*(x,r) cp,*(y,z), / = 1-3 , 

*■' (12.0,*) 

6 

<Pf/ = *» «=l-3» y -1—6 . 

*=1 


This displacement field will be used as the basis for the remaining development, where no limiting 
assumptions are made corresponding to load type. 


Ordering Scheme 

To identify and delete higher order terms which are produced during the derivation of the governing 
equations, it is necessary to define an ordering scheme. First of all, let assume that the beam is 
relatively long and slender such that the geometric ratios of the cross section to the length (y/L, z/L) are 
in the order of e which is defined to be equal to 0.2. Since x is the centroidal axis of die beam, {xJL , 
L d( )/dx ) are in the order of one. To assign the orders of others terms, a study of the deflection 
patterns of a tip loaded anisotropic cantilever beam composed of material properties of High-Strength 
Graphite/Epoxy (see Table 1) has been performed using the exact results obtained by Ie and Kosmatka 

(1992.b). Due to the fact that the beam spins, the order of uIL has been increased up to e 2 . 


Table 1: 

Typical Material Properties of High 

*11 

145.0 GPa 

II 

cs 

10.0 GPa 

G 12 = G 13 

4.8 GPa 

v 12 “ v 13 

0.25 

v 23 

0.40 

P 

1580 kg/m' 


Taking only up to the third order of s, an ordering scheme will be obtained in the following manner 
Order 1 : 

!• L V- < 13 *> 

Order e: 


6.5 


Z z. £ 

L ’ L ’ L ’ 


e,. e y , e; 


(13.6) 


Ordered 


Order e3 : 


it 

L’ 


d(¥l. M>2, VP3) 

ay 


(Vl» V2, H>3 ! ) 


0(vi» M>2> V 3 ) 
dz 


d(Vu V3) 

dx 


(13 jc) 


(13 .d) 


Nonlinear Strain-Displacement Relations " 

A geometrically nonlinear strain-displacement relationship is being adopted as follows 


dx 2. 

/du\ 2 + /dvi\ 2 fdw\ 2 
\dx) \dx‘f \dx] 

1 

£yy - 

E =^+I 

Ezz dz 2 

( + /dv\ 2 / dvv j 2 
\dzj \dzf \ dz / 

> 

y>z = 

_ dujdw + du dn + dv dv + dw dvv 
^ zx dz dx dz dx dz dx dz dx 

Yxy = 


dv 

dy 


4l-F4 

21W 


dv\ 2 /dvv' 2 

M W 


dvv + dv + du du + dv dv + dvv dvv 
dy d z dy dz dy dz dy dz 

dv + d« + dw du + dv dv + dvv dvv 
djc dy dx dy dx dy dx dy 


(14.*/) 


Substituting eqns (l.a-c) into the above equations and applying die ordering scheme as shown in eqns 
(13 .a-d), the strain array can written as follows 


{eHXlUl , 


(15) 


where the matrix [A] only depends upon cross sectional dependent functions and { € } is an array of 
macroscopic strain measures which depends on the six variables u>v,w , 0*, Qy and 0^ . [X] and 
{ € } can be decomposed into linear and nonlinear parts as follows 


[xri [x,] [x 2 J [x„ 3 J ] , 
(e }={ € /}+ { €n/} + (eL). 

See the Appendix for detailed expressions of [X] and { € } . 
Strain Energy 

The strain energy is defined as follows 


(16) 

(17) 


"■if i ,f 


■} T [C\{t)dA dx . 

Substituting eqn (151 into eqn (18), the strain energy can be written as follows 


( 18 ) 


6 . £> 



(19) 


. . U = l( 

where the matrix [Se] is called section constant and is defined as follows 

[Se}= j [X] T [C] [X)dA . 

Furthermore, substituting eqn (17) into eqn (19), the following will be obtained 

U=\ f (! € ^) r+ { € «/) r+ ( € /) 7 )[ 5 c l(( € ri + ( € ^i + { € »»)) <ix 

1 Jo 


(20) 


(21) 


Finally, cany out the multiplication, neglecting the term associated with multiplication between the 
third order strain measures, the following expression for the strain energy will be obtained 


U = 2 j ( ^ € ( € /} + 2 { € /} 7 [5e] { € li\+ 2 { € /) r [Se] { € ^/} 
+ { € [5e]{€^)+2 [5e]{e^/})cix . 

Kinetic Energy 

The kinetic energy of the beam is defined as follows 


< I. 


p V.V dA dx. 


where p is the mass density and V is the velocity vector which is equal to 




in which r is the position vector defined as 


r- {H x , H y , H z ) j j y j+U+u, y+v, z+w)[ e 

ul 

Transforming into the beam coordinate system, the position vector can be written as follows 


r= (h x +x+u, hy+y+v , A z +z+vvV 


c y / » 


( 22 ) 


(23) 


(24) 


(25) 


(26) 


6 . 7 


where (h x ,h y ,hj = (H x ,H y ,H£ [B]T, and [B] is a transformation matrix between the two sets 
of the coordinate systems by rotating through Euler angles Pj, 02> P/ which [B] = [B^] [B 2 ] [B 3 ]. 


Bi=4 


1 0 0 
0 cos Pi sin Pi 
0 -sin Pi cos Pi 


b 2 4 


( cos P 2 0 -sin P 2 
0 1 0 


b 3 H 


cos P3 sin P3 0 

-sin P3 cos P3 0 
0 0 1 / 


(27 .a-c) 


\ sin P 2 0 cos P 2 / 

Substituting eqn (26) into eqn (24) and making use of eqns (27.a-c), the velocity vector becomes 
V= (u+Oy^+z+wJ-a^hy+y+v))^ + (v-n x {h z +z+w)+n£h x +x+uj)e y + 

(w+Cl x (h y +y+ty-Cly(h x +x+u))e z , 


(28) 


in which ( Cl x , Cl y , Q^) = (0, 0, Q ) [B ]T . 


The Governing Equations (Equ ations of Motion) 

The six governing equations for the beam are determined by applying the Hamilton's principle as 
follows 


6 



(-U+T+W e )dt=0 




where V, 7*and W e are the strain energy, kinetic energy and work done by external forces, 
respectively. 

The variation of the strain energy can be written as follows 



dx. 


(29) 


(30) 


Substituting the expression for the strain array from eqns (15) and (17) into the above equation, the 
following result will be obtained 



where {B} is an array of stress resultants with forty eight elements defined as follows 


{*) 


I. 


[Xf{o)dA. 


(31) 


(32) 


The variation of the kinetic energy as shown in eqn (23) can be decomposed into three parts as follows 



(Br 2 +6r 1 +6r 0 )</r, 


(33) 


6.8 


where 



(34 .a-c) 



The above equations are associated with inertia, Coriolis and centrifugal forces, respectively. Note that 
{ * } indicates the derivatives with respect to time. The array {< 7 } is defined as follows 

{?} r = («, v, w, e x , e y , e z , u, v\ w’, e„ 0 y , e z ) . (35) 

Matrix [m\, [m c ] and are the mass, Coriolis, and centrifugal softening -type matrices, 
respectively. Arrays {f co } and {f c f} are the Coriolis and centrifugal force-type arrays, respectively. 

Due to the fact that Q is constant, the array {f co } will not contribute to the governing equations. The 
details of the expressions can be found in the Appendix. Integrating by parts the expressions in eqns 
(34. a-c) with respect to time, the following results will be obtained 



(36. .a-c) 


where 

{Z)=[m]{q}, {Z\ = [m co ]{q), (z) = -[* s ] [q ) . (37.o-c) 


in which [m co ] = -2 [m c0 \. 

The variation of the external work done is as follows 


6W,= 



(Px + Py 6 v + p z 5w) dx . 


(38) 


Substituting eqn (31), eqns (3 6.a-c) and eqn (38) into eqn (29) and taking any necessary integration by 
pails, the boundary conditions and the six governing equations will be obtained as follows 


Boundary Conditions : 

Specify u or N u which is equal to 




N u 9 i?i-i? 7 + /?i 9 V , +/f 25 H’'+i? 3 i 9 I + i? 370 y-^ 430 z + ^ 7 + ^ 7 + 27-/7 • 
Specify v or N v which is equal to 

N v *R2-R'i+Rl9U'HRl3+2R20)vHR21+R26)W ,+ (Rll + R32)Qx + (R2l + R3s)Qy- 

(i?2O + ^44)®z + ^220* + ^23®y"* - '^24®z + ^8 + '^ + ^8'^ • 

Specify w or N w which is equal to 

N w »/?3-/?9+/?25K' + (/?21 + fl26) v+ (^13 + 2f?27)w' + (^18 + ^33)8;c + (^27 + ^39)0)r 

(f?26 + ^45)0z + ^28®* + ^290y + ^3()Qz + ^ + ^ + ^'^ • 


(39 .a) 


(39 .b) 


(39 .c) 


Specify ©x or M x which is equal to 

M x ■ fVfiio + ^22 v + fl28w' + /?340* + tf4oOy-fl460^ + Zlo + Zio + 2lo-/io• 
Specify 0 y or My which is equal to 

M y ■ J l ?5+f?9-^i 1 +f?23V ,+ i?29W'+f?350x + i l ?410y-f?47©z + ^ll + ^ir , ' Z ll-/n ■ 
Specify 0 Z orM z which is equal to 

M z ■ f? 6 -^ 8 -^i 2 + ^ 24 V ,+ /? 3 OW'+i? 360 x + ^ 420 y-^ 480 z + Zl 2 + Zl 2 + Zl 2 -/l 2 • 


(39.40 


( 39 .e) 


( 39 ./) 


Specify u\v\w\Q 0 y and 0 ' z or M u , M v , M w , R x ,*y , and > respectively , which are 
defined as follows 


( 40 .O-/) 


M u m Ri , M v ■ /?g > M w -/? 9 , 

#**^10, fly"/?ll, /?z"^12- 

The six governing equations : 

^--(Zi+Zl+Zr/il+p^, ^-(Z2+^+Z 2 -/2) + /V=°» ^-(Z3 + Z 3+ Z3-/ 3 )^z=0, 

(41.0-/) 

N v -lZ6+%+%-/d=0, 


^UJV-(Z4+Z4+Z4-/ 4 M) , ^-F.-fZs+Zs+Zs/sM), ajc 


5M Z TT 


where 


F- -f?3iU'-(/?i7 + ^32)v'-(^18 + ^3^'-(^14 + ^15)0x-^330y + ^32Oz + ^35Oy-^360z , 

F v ■ /?2 + K43« Mtf2O + K44)vMfi26 + rt45)w' + i?32®* + (rtl6 + tf38 + K45)0y- 

(/?14 + 2/?44)0 z + ./?460x + f^47®y > 

F w ■ f? 3 +i ?370 + (i?21 + f?38) v + (^27 + ^39) >v + f?330x + (^lS + 2/?39)0y _ 

(f?16+i?3g+/i45)0Z + i l ?4O0* + /f42Oz * 


( 42 . o-c) 


£>. 10 
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S ELEMENT MODEL 


The finite element model is developed using the weak form formulation. The six variables u, v,w 
fix. 9y and 0 Z are interpolated using the expressions of the form 


4 5 

U = £ Uj <fc" , V = £ V, <fc V , 

r =1 i=l 


W 


= 1 w,4*\ 

1=1 

6, = £ e'*', e,= £ e/V, e ; =£ 


1=1 


4 

i 

i =1 


(43 .a-f) 


»=i 


where the variables ufix, 6y and 0 Z are interpolated using four-node Hermite polynomials and v and 
w are interpolated the same way as u except that one additional node being added in the mid-length of 
the elements. Substituting the above equations into eqn (29), the associated governing equations (the 
equations of motion) related to the nodal displacement {Q} can be written as follows 

+ . ( 44 ) 

where 


{Q) T = { (Ml, « 3 , W 4 ), (Vf, v 2 , v 3 , V 4 , V 5 \ (Wi, W 2 , w 3 , W 4 , W 5 ), 

(e, x . el el e 4 *), (ef, el el ei'). (ef, e|, el e 4 ') } , 


(45) 


[M] is the mass matrix, [C] the gyroscopic matrix, [AT/ ] the linear stiffness matrix, [Krij ] the nonlinear 
stiffness matrix, [K c f] the centrifugal stiffness matrix, {F c } and {F^ the force vectors related to the 
external and centrifugal loads, respectively. 


NUMERICAL EXAMPLES 


In this section, numerical examples are presented for beams with elliptical cross section composed 
of an anisotropic homogeneous material (see Fig. 3.a) where (a) and (b) are one-half of the major and 
minor dimensions of the cross section, respectively. Through out these examples, a value of (b/a) 
equal to one half has been chosen. Furthermore, the beam is assumed to be composed of a single set of 
unidirectional high strength graphite/epoxy fibers, which is a transversely isotropic material that has 
five independent constants (see Table 1) 



Fig. 3. a A tip-loaded cantilever beam 


To achieve orthotropic and/or anisotropic beam behavior, the material reference frame (1,23) is 
oriented relative to the beam Cartesian coordinate frame (x,yj) using reference angles (P) and (a), 
which are defined as rotations about the positive x-axis and the positive material 3-axis, respectively. 
See Fig. 3.b, where the transformation relation between (133) and (x,y f z) is equal to 


IM 

cos(a) 

sin(a) cos( P) sin(a) sin(fi) 

X 

2 - 

- sin(a) 

cos(a) cos( P) cos( a) sm(P) { 

y 

1 3 / : 

0 

- sin( P) cos(P) 

,Z , 


The resulting 21 unique material compliance coefficients (Sjj, i j=l-6) are determined using standard 
transformation techniques (Lekhnitskii, 1963). Numerical examples for static and free vibration cases 
are presented for the case of a = P = 30 degrees, L/(2a) = 10. The model predicts the coupled 
displacements and stresses exactly. 

Example 1 : Static. Nonspinning 

Consider a cantilever beam fixed at the root with a single tip load P z (see Fig. 3. a). Unlike the 
isotropic counterpart, results show that all six stress components are present Fig. 4.a-j show the 
countour plots of these stresses. 

Normal Stress : 

Fig. 4.a-d show the changes of distributions in the cross sections at different values of x. 
Dotted lines show negative stresses. At the tip (Fig. 4. a), unlike the isotropic counterpart (in which 
this stress is equal to zero), ct** distributes in a such a way such that the resultant is equal to zero. 
This stress will gradually becomes less dominant as the value of (L-x) increases (for cross sections 
closer to the root). Fig.4.b shows the stress distribution for a value of {L-x)~ b, Fig.4.c for (L-x)= 

2 b, and Fig.4.d for (L-x)=L (at the root) where the stress distribution is dominated by the "classical’' 
normal stress distribution. 

In-Plane Stresses &yy » oy? and t 

Due to the anisotropy of the material, these stresses exist Because the shear forces through out the 
length are constant these stress distributions do not change with respect to x . As they have to be, 

Oyy and (see Fig. 4.e,f, respectively) vanish at two locations, i.e. at y = ±a , z= 0 for Oyy and 
at y = 0 , z= ±b for . Likewise, Oy Z vanishes at four locations, i.e. at y = ±a , z= 0 and at y 
= 0 , z= ±b (see Fig.4.g). 

Transverse Shear Stresses && and o^-y : 

Fig. 4.h,i show the countour plots of transverse shear stresses and oj y , respectively. Fig.4.j 
shows a qualitative vectorial scaled plot of the resultants between these two stresses. From this figure, 
it is interesting to see that the stress distribution clearly does not resemble that of an isotropic 
counterpart (i.e. a paraboloid with a maximum occurring at the centroid). Instead, the distribution is 
nonsymmetric with a maximum occurring on the outer edge. 



x 



Fig. 4 Stress distributions of a tip-loaded cantilever beam. 
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Example 2 : Free Vibration. Nonspinning 

For free vibration analysis, consider the same beam as in example 1 but with free-free boundary 
condition. Table 2 shows the mode number, the values of Xj (no dimension) and information about 
the mode shapes, where 


ft 


x< .p ^Un. 

2 7t L 2 V p A ’ 


(47) 


fi is the i-th frequencies (in hertz), and lyy is the moment of inertia with respect to t hey ax is. As 
expected, the first six modes represent rigid body modes. For lower frequencies, couplings exist 
mainly between bending and torsion. Complete couplings (bending, torsion, shear, and extension) 
exist for higher frequencies. 


Table 2: 


Mode 

Number 


1 to6 

0 

7 

8.33 

8 

16.42 

9 

22.66 

10 

42.93 

11 

43.40 

12 

49.99 

13 

69.80 

14 

80.69 

15 

94.53 

16 

99.45 

17 

101.58 

18 

125.76 


Information about the mode shapes 

Rigid body motions 

First bending in the x-o-z—plane with little torsion 
First bending in the x-o-v— plane with l ittle torsion 
Second balding in the x-o-z-plane with little torsion 
Third bending in the x-o-z—plane with some torsion 
Second bending in the x-o-y— plane with some torsion 
First torsion with some third bending in thex-o-z-plane 
Fourth bending in the x-o-z— plane with some torsion and li ttle extension 
Third bending in the x-o-y-plane with some torsion 
First extension with some torsion and little bending in the x-o-y— plane 
Second torsion with little bendings in both planes 
Fifth bending in the x-o-z- -plane with some torsion 
Fourth bending in the jr-oy— plane with some torsion, little extension and 
little fourth bending in thex-o-z-plane 


CONCLUSION 

A formulation of a nonlinear theory for spinning anisotropic beams having arbitrary cross-sections 
has been developed. The di splacement field is assumed to compose two parts, i.e. the standard 
kinematic relations (expanded for a three-dimensional case) to describe the global beam behavior and a 
supplementary additional field that represents the local warpings within jmd out-of the cross-section 
plane. It has been shown that, in the most general case, these beam warping functions may be as 
many as eighteen. Furthermore, it was assumed that the magnitude of this additional field is directly 
proportional to the local six stress resultants. 

Using a weak form finite element based numerical technique, preliminary nu meri cal results for a 
static and a free vibration cases have been obtained. These numerical examples show that the model 
can predict the displacements and the stresses exactly for a shear tip loaded cantilever beam. Due to the 
anisotropy of the material, stress distributions do not resemble the isotropic counterpart for all six 
stress components. Numerical results on the free vibration case (a free-free beam) show the coupled 
mode shapes of bending, torsion, shear and extension, especially for higher frequencies. 
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APPENDIX 

The matrix [X ] is decomposed into three parts as follows 


1 <P21,2<P31,3 

<P21,3 + <P31,2 

<Pll,3 

<Pl 1,2 

- <P22,2<P32,3 

<P22,3 + <P32,2 

<Pl2,3 

1+<P12,2 

- <P23,2<P33,3 

<P23,3 + <P33,2 

1+<P13,3 

<Pl3,2 

- <P24,2 <P34,3 

<P24,3 + <f>34 ) 2 

y + <Pl4,3-Z + <Pl4,2 

Z <P25,2«P35,3 

<P25,3 + <P35,2 

<Pl5,3 

<Pl5,2 

-y <P26,2<P36,3 

<P26,3 + <P36,2 

<Pl6,3 

<Pl6.2 

<Pll - 

- 

<P31 

<P21 

<Pl2 - 

- 

<{>32 

<P22 

<Pl3 “ 

— 

<t>33 

<P23 

<Pl4 “ “ 

- 

<P34 

<{>24 

cpis - 

- 

<P35 

<P25 

<Pl6 ~ 

- 

<P36 

<P26 _ 


1 


. [X. 2 /] 
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[XnM r = 


— — 

— 

— 

921,3 921,2 

- - 

- 

- 

922,3 922,2 

- - 

- 

- 

923,3 923,2 

-z - 

- 

- 

924,3 924,2 

- - 

- 

- 

925,3 925,2 

- - 

- 

- 

926,3 926,2 

- - 

- 

- 

931,3 931,2 

- - 

- 

- 

932,3 932,2 

- - 

- 

- 

933,3 933,2 

y - 

- 

- 

934,3 934,2 

- - 

- 

- 

935,3 935,2 

- - 

- 

- 

936,3 936,2 

- 9*31,2 

-<P21,3 

9*31,3“9*21,2 

— 

— 

- 9*32,2 

"9*22,3 

932,3-922.2 “ 

— 

- 9*33,2 

"9*23,3 

933,3-923,2 “ 

— 

“ 9*34,2 

“9*24,3 

934,3-924,2 Z 


“ 9*35,2 

“9*25,3 

935,3-925.2 - 

— 

" 9*36,2 

"9*26,3 

936,3"926,2 - 

— 

- - 

9*11,3 

9l 1,2 

1 

— 

- - 

9*12,3 

9l2,2 

— 

— 

— - 

9*13,3 

9l3,2 

— 

— 

- - 

9*14,3 

9l4,2 

— 

— 

- - 

9*15,3 

9l5,2 

z 

— 

- - 

9*16,3 

9l6,2 

-y 

— 

“ 9*11,2 

— 

9l 1,3 

— 

1 

- CPl 2,2 

— 

9l2,3 

— 

— 

- 9*13,2 

- 

9l3,3 

— 

— 

“ Vl 4,2 

- 

9l4,3 

— 

— 

“ 9*15,2 

- 

9l5,3 

— 

Z 

. - <Pie,2 

- 

9l6,3 

— 

-y _ 


respectively. Similarly, the array { € } is also decomposed into three parts as follows 

{ € lY"— (u , V - 8 *, W + 0 y, 8 x , 0 y, 0 z , u , V “ 0 Z , W + 0 y» 0 JCj®y ^ 

0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) , 

{ e J ; ) r = (0, 0, 0,0,0, 0, 0, 0,0,0, 0, 0, 0.5 [(vf + (wf], 0.5 [(0 x f+ (e #1 

0.5 [(0 x f + (0y.Fl -0y Qz, V' 0 X , w' 0 X , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 

0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) , 

( € lij T = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, v V, vtv'-0 z ), 
vtw'+Oy), v'0 x , v’0y, v'0 z ,wV, wtv'-ej, wtw'+0y), 
w'0 x , w'0y, w'0 z , 0 x u', 0 x (v'-0 z ), 6 x (w'+0 y ), 0 X 0 X , 

0 x 0y, 0 x 0 z ,0yM , 0y(v ~0 Z ), 6y(w + 0yl), 0y6x> ®y®y* ®y®z» 

-Q z u\ -0*(v -0j, -O^W + 8y), -0 Z 0x, “Oz©y* *^z®z) • 

Matrices [m] and [m c ] are defined as follows 

6. It 


(A. lx) 


0 , 

(A. 2. a) 
(A. 2. b) 

(A.2.c) 


where 


["*]=j f P ({mi } {mi } r+ {m 2 } {m 2 } 7 +<m3} {m3} 7 ) dA , 

["» c ]=f P {O x (i{m2}{m3) 7 '-{m 3 }{m2) 7 )+Qy({m3}{mi} r -{mi}{m 3 }^ 

i){w 2 ) r -<m 2 }{iiii} 2 )} </A, 


{mi) r =<l,0,0,0,z+cpi3, -y-<Pi2, <Pi 1, <Pi 2 , «Pi 3 , <Pi4, <Pi5, <Pie) , 
{m 2 } 7 ={0, 1 , 0, -Z, <P23, -Cp22, <P21, <P22, <P23, <P24, <P25, <P2c) , 
{m 3 } r =(0,0, l,y, CP33, -<p32* <P31> <P32» $>33> <P34> <P35> <P36) • 


(A. 3. a) 


(A.3.b) 


(A.4.a-c) 


Matrix [fc 5 ] is defined as follows 


where 




(A. 5. a) 


&) = j p { (Oy +0/) {mi) {mi } r+ (^ + ^ 2 ) ( m 2 ) {m 2 } r+ (n^+Q^) {m 3 } {m 3 } 7 

-20yQ z {m2} {m 3 } r -20 x 0.{mi} {m 3 } 7, -20 I 0 y {/7ii} {m 2 ) r } dA . 

The array {f c j$ is defined as follows 

{fcjhj p { (n y 2 +nf)(/ lj +jc){m 1 5+(n J ?+ni i )(/. y +^){m 2 i+(oi i +a y 2 )(/ l2 +z){m3}- 

ft*(^Vz)+Oy(Vy)) < m i ) A'( n *( V*) + ft*(A* + *)) { m 2 )- 

toz(tox(hx + x)+Cly(h y +y)){m 3 ) } dA . 


(A.5.b) 


(A. 6) 
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Abstract 

An analytical model is developed for assessing the 
extension-bend-twist coupling behavior of nonho- 
mogeneous anisotropic beams with initial twist. The 
model is formulated as a coupled two-dimensional 
boundary value problem, where the displacement 
solutions are defined with pretwist-dependent func- 
tions that represent the extension, bending, and 
torsion, and unknown functions that represent local 
in-plane deformations and out-of-plane cross-section 
warping. The unknown deformation functions are 
determined by applying the principle of minimum po- 
tential energy to a discretized representation of the 
cross section. Numerical results are presented that 
fully verify this approach and illustrate the strong 
extension-twist coupling behavior present in 
pretwisted beams with thin-wall laminated composite 
cross sections as a function ol ply angle, initial twist 
level, and initial twist axis location . Cross-sections 
analyzed include; thin laminated rectangles with ei- 
ther asymmetric or symmetric ply stacking se- 
quences and a thin-wall single cell D-section com- 
posed of a graphite/epoxy woven cloth. 


introduction 

From tilt-rotor aircraft to jet turbines, rotor 
blade manufacturers are incorporating fibrous com- 
posite materials into their current designs as a 
means of reducing weight and costs, and controlling 
deformations. In a general sense, a laminated com- 
posite rotor blade can be described as an elastic 
beam that exhibits generally anisotropic behavior, 
where it's shape is generated by rotating a nonho- 
mogeneous airfoil (irregular) section about an initial 
twist axis. The beam can have a helical line of cen- 
troids, since the section centroid is not required to 
lie on the initial twist axis. Thus, the application of a 
simple extension load will result in elongation, 
bending, and twisting of the beam. This coupled 
behavior is dependent not only upon the material 
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property definition of the laminated section, but also 
upon the initial twist axis location and the initial 
twist level. 

Exact analytical solutions to this type of 
problem do not exist and are generally intractable. 
Fundamental studies [1,2] derived the governing 
two-dimensional coupled equations for the exten- 
sion-torsion behavior of pretwisted and spiral 
(helical) isotropic bars. Closed-form solutions for 
the extension-torsion behavior of helical isotropic 
bars with simple (off-centered circle and ellipse) 
cross sections were developed in [3] using a dis- 
placement formulation. Numerical results clearly illus- 
trate the interaction between pretwist and local 
cross-section deformations. A recent study [4] de- 
veloped an analytical model for pretwisted isotropic 
beams with an arbitrary cross section, where the 
Ritz method is applied to determine the pretwist de- 
pendent in-plane deformations and out-of-plane 
warping of the cross-section before studying the 
extension-bend-twist coupling behavior. Numerical 
results demonstrated the pronounced effects that 
pretwist and initial twist axis location have on the 
section deformations, extension-torsion coupling 
behavior, and section properties of solid and thin- 
wall multi-cell airfoil sections. Investigators have 
developed other isotropic models based upon either 
thin shell theory [5] or approximate technical beam 
theories (for example, [6-10]). 

Pure bending of pretwisted isotropic bars with 
simple homogeneous sections has also been 
addressed by investigators. Maunder and Reissner 
[11] developed approximate solutions using a thin 
shell theory for narrow rectangular cross sections. 
Goodier and Griffin [12], using a stress formulation, 
developed an elasticity model assuming that the so- 
lution can be represented by a pretwist dependent 
power series. Results using the first few terms of 
the series for a thin elliptical cross section show 
that curvature will increase significantly for the stiff 
plane of the cross section, but the curvature in the 
soft plane remains virtually unchanged. 

Independent of the above research on 
pretwisted isotropic beams, investigators have de- 
veloped solutions for the behavior of prismatic 
anisotropic beams with a nonhomogeneous irregular 
cross section. Initially, a mathematical formulation 


Tl 


with an existence proof was derived based upon an 
assumed displacement field, but no numerical results 
were given [13J. Approximate solutions, which in- 
volve solving a coupled two-dimensional elasticity 
model via the Ritz method, have been developed 
[14-16]. In [14], the solutions are determined by un- 
coupling the local cross-section deformations from 
the global beam deformations and solving both 
simultaneously, whereas in [15, 16], the global beam 
solutions are derived first using Saint-Venant's 
inverse method and then only the local in-plane and 
out-of-plane section deformations require 
calculation. 

The objective of this paper is to develop an 
analytical model for studying the extension-bend- 
twist coupling behavior of nonhomogeneous 
anisotropic beams with initial twist. The model is 
formulated as a coupled two-dimensional boundary 
value problem, where the displacement solutions are 
defined in their most general form including: (1) 
pretwist-dependent functions that represent the ex- 
tension, bending, and torsion, and (2) unknown 
functions that represent the local in-plane deforma- 
tions and out-of-plane warping of the cross section. 
The unknown deformation functions, which are as- 
sumed to be proportional to the local axial strain, 
bending curvatures, and torsion twist rate, are de- 
termined by applying the principle of minimum po- 
tential energy to a discretized representation of the 
cross section (Ritz method). Finally, the extension- 
bend-twist coupling behavior is studied using the 
equilibrium equations of the cross section. This 
model has direct applications to both highly 
pretwisted aviation propellers and jet turbine (turbo- 
fan) blades, which have thin built-up solid laminate 
sections, and composite tilt-rotor blades, which 
have thin-wall closed-cell laminate sections. 


Consider a long elastic beam, of length I, where 
the lateral surface is generated by rotating an 
arbitrary nonhomogeneous cross section about an 
initial twist axis (z-axis). This cross section is de- 
fined using (n) triangular and/or quadrilateral subre- 
gions, where each subregion can have a unique ho- 
mogeneous anisotropic material definition. See Fig. 
1. The beam may have a helical line of centroids 
since the modulus weighted section centroid is not 
required to lie on the initial twist axis. A space-fixed 
orthonormal vector set ( x,y,z ) and a curvilinear co- 
ordinate system (£,tj ,z) are used to analyze the 
beam, where £ and rj align with the x and y axes at 
the beam root (z=0), but rotate with the section 
about the initial twist axis. Thus, both the section 
geometry and the material properties are functions 
of | and r\ only. The two coordinate systems are 
related using 

I g \ cos(az) sin(az) 

\ ij | -sin(az) cos(az) 



where a is the initial twist per unit length. The con- 
stitutive relations for the i th subregion of the cross 
section, are given as 



C 


4e»{ . 

J l 1 




(2 .a-b) 


where the i th stress and strain vectors are given as, 


Three sets of numerical results are presented. 
Initially, the extension-twist behavior of a flat 
(untwisted) laminated plate with an asymmetric 
(angle ply) lay-up is analyzed to verify that the cur- 
rent approach reduces to classical laminated plate 
theory. Second, the extension-bend-twist behavior 
of a thin solid laminated strip with initial twist is 
studied for different ply angle orientations and 
stacking sequences (asymmetric, symmetric), 
pretwist levels, and initial twist axis locations. 
Finally, a pretwisted beam having a thin-wall D-sec- 
tion composed of a graphite/epoxy woven cloth, 
similar to the effective structural section of a tilt-ro- 
tor blade, is analyzed to illustrate how ply orienta- 
tion and initial twist can be combined to produce ei- 
ther maximum or minimum extension-twist coupling 
behavior. 


Theoretical Derivation 


M- 


(2.C) 


| ®zz> ^ 7Z^> T qz^’ 


w- 


(2 .CT) 


{ Hi {i) ’ C W' } • f ^ (,) ’ y nz 0) ’ X-z (,) • 


and the material stiffness [CW] and compliance 
[SO)] matrices for each of the subregions in the 
curvilinear frame (<-,jj,z) must obey [C^] =[S^]‘ 1 . 
These matrices will be fully populated with up to 21 
distinct coefficients when the subregion material 
classification is either anisotropic or the subregion is 
composed of fiberous composite materials, where 


the principle fiber directions do not align with any of (,) _ _ nfl w (;) 

the curvilinear coordinates. 01 ' oV ' ^ 


The applied stress distribution on the ends of 
the beam (z=0,f) is statically equivalent to an applied 
extension force P that acts along the initial twist 
axis arid a general moment M that can be de- 
composed into components M x , Mu and M z about 
the x, y, and z axes, respectively. Furthermore, the 
ends of the beam are free to warp so that the twist 
is uniform along the length. A general moment is 
required for two reasons; (1) P does not act 
through the centroid and thus effective bending 
moments are produced, and (2) the solutions to the 
bending and torsion problems are coupled for gen- 
erally anisotropic beams due to the presence of the 
C34 and C35 terms in the material stiffness matrix. 
Assuming that the body forces are negligible and a 
stress free condition exists along the lateral surface, 
then the stresses within the cross section must 
satisfy the following equations of equilibrium; 


I 

/-I 


I OzzM dA M = P , 



v(') » v 0 (z) +q9 0 (z) + \ff 2 W . (5.a-c) 

wM = W 0 (Z) -Z0T) (Z) + T 10§(2) + V 3 W , 

where uo, v 0 , w 0 represent z-dependent displace- 
ments in the tj, and z directions, respectively, 4>£, 
<t>q, and 9o are rotations of the cross section plane 
about the |, tj. and z axes, respectively, \y\ and W2 
are deformations in the section plane (including 
Poisson contractions), and y/3 describes warping 

out of the section plane. These functions 
V2^ l \¥3^) are assumed to; (1) be directly propor- 
tional to the axial strain, bending curvatures, and 
twist rate within the cross section, (2) be uniform in 
the curvilinear coordinate frame (function of | and 77 
only), and (3) have first-order continuity across the 
subregion boundaries. 

Assuming a two-dimensional strain state, one 
can derive the final form of the z-dependent 
functions as 


X f T]<?zz {i) dA M = Mz , 

/-I V 

(3 .a-d) 

X [ dA ( ' )= - M T) . 

'-I V 


u 0 {z) = — [(az)cos(az) - sin(az) j 
a 2 1 


K C t I 

— 1 * cos(az) - (az)sin(az) 

2 * 1 


a 


X - »jr* 2 W| dA W = M z , 


Mz) = 



(az)cos(az) - sin(az) 


I 

I 


where, A^) is the area of the i th subregion and M* 
and Wtj are components of the applied moment 
about the ^ and 77 axes that satisfy 



r 


- 


cos (az) 
-sin(az) 


sin(az) | M x \ (4) 

cos (az) \ M y j 


The displacement distribution for each subregion 
can be written in its most general form as global 
functions that represent extending, bending, and 
twisting of the beam and local functions that 
represent warping of the nonhomogeneous cross 
section; 


+ — jl - cos(az) - (az)sin(az)j , 

,2' ' 

( 6.a-/) 


a‘ 


w Q {z) = e z , 


<Mz)=— jl - cos(az)) + ^jsin(az)) , 
a | j a \ j 

O n( z) = -^ji - cos(az)j + -^jsin(az)J , 


9 0 (z) « 9z , 


where, e, k~, k tj, and 9 represent the extension 
strain, the bending curvatures of the beam in the |-z 
and tj -z planes, and the elastic twist per unit length, 
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respectively. The current approach reduces to the 
models in (3, 4, 12] for isotropic materials. 

The i th strain components are found using the 
displacement relations from (5.a-c) and (6. a-f) 

= e - + riK-q + , 

7 ^ = 0 * + +<#¥'■ i^+Oi^W), (7. a-/) 

= - 0 | + 1/3.^ - a( y/^ -OV'i (/) ) - 

W' )= % ( ' )+ ^ 2 .= ( ' ) * 

where, the symbol D is an operator defined as 


d=t)S--=£- 

51 sn 


( 7 . 9 ) 


Determination of the Local Cross Section 
Deformations 

The local deformation functions for an arbitrary 
nonhomogeneous anisotropic cross-section are 
determined based upon the principle of minimum 
potential energy along with a discretized represen- 
tation (finite element modeling) of the cross section. 
Although the displacement field is fully three-dimen- 
sional (Eq. 5.a-c), it is explicit in the z direction, thus 
only the two-dimensional cross section needs to 
analyzed. This approach has been applied success- 
fully to such problems as torsion and flexure of 
prismatic isotropic [17], monoclinic [18], and 
anisotropic beams [1 5, 1 6] , and the extension-tor- 
sion coupling behavior of pretwisted isotropic 
beams [4], 

The local deformations must be determined for 
each of the four cases (viz. extension, bending 
curvatures in the £-z and 77-z planes, and the elastic 
twist rate). Standard isoparametric finite element 
methodology is employed so that most of the details 
can be omitted. Each material of the cross section 
is approximated using either quadrilateral or triangu- 
lar subregions where the local deformations are rep- 
resented as 


V'i 


(') = 


V2 


(/) = 


Vs 


(') 




N^d 77) 


N^dn) 


\ 1 r 



(8 .a-e) 


L«\ 
J\ 3 1 ' 


where [/v^(!,r ?)] is a bi-quadratic isoparametric in- 
terpolation function and { */ / i (0), and {^3^} 

are nodal displacements on the i}" subregion 
boundary in the 77, and z directions, respectively. 
The strain vector (Eq. 2 .d) of the ft subregion can 
be written in matrix form in terms of the unknown 
local deformations and the extension strain, bend- 
ing curvatures, and elastic twist rate by substituting 
the interpolation functions (Eqns. 8.a-c) into (Eqns. 
7 .a-/): 


j f (')J = S b W 


( \ , r 1/ \ 


m 

\ j [ '■ri 

(9) 



gg 


where 


N li) dv)-; 


B 


U) 


aN U) {*,T}) 

aDN l \q , n ) 

/Ac,??).* 


N lh dri). n 


aDN'' ) d rj) 
-aN'^dTj) 


aDN l '\-,n) 


(10. a) 


jtP w J = 1 1 j , { ^ j ’ | J J (10.fr) 


and 


0 0 0 0 

I 0 0 0 0 

i •; i o 

o o o 4 

I 0 0 0 -)] 

j_0 0 0 0 


(10. c) 
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jb = !©. fj. Kjj, 6 j 


(10. d) 


Similarly, the displacements (Eqns. 5 .a-c) could also 
be written in matrix form in terms of the local cross 
section deformations and the vector { b }. 

The principle of minimum potential energy is given as 

sn = X 5u{!) * 5v ^e = 0 (11) 

ti 

where n is the number of subregions, 8 u( 0 is the 
variation of the strain energy with respect to the 
unknown local deformations of the ft subregion 
given by 


dA^dz, 


(12. a) 

and 5Wg(') is the variation of the work of external 
forces of the ft subregion that results from the 
applied tractions on the beam ends; 
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and the force matrix is presented as 
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Since both the stiffness matrix ( [K(')\ ) and the 
force matrix ( [F^O)] ) are linearly dependent upon 
the beam length (L), then the calculated local de- 
formations functions are length independent and 
( L ) can be dropped from the above equations. 

The matrix equations of Eq. (13) are assembled 
into a complete model of the cross section using 
standard finite element procedures. Unit solutions 
for the local deformations (V'i ,¥ 2 ,^ 3 ) can be 
calculated for each of the four cases of { b ) by 
setting the appropriate value in the array {b} equal to 
unity and the remaining three to zero. Thus, the cal- 
culated deformation functions can be written in ma- 





trix form as 


SWe - 

| (<1, ('1 «V- !'1 </). ('1 1 
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+c zz d V3 
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(*• 0 ) 


dA 1 * - 


The virtual work expression will reduce to zero since 
both the stresses and the local cross section de- 
formations are assumed to be independent of the 
axial coordinate (z). A set of linear algebraic equa- 
tions for determining the local cross section defor- 
mations in terms of { b } is obtained by substituting 
Eqns. (9) and (12. a) into Eq. (11) and carrying out 
the integration over the beam volume. Writing this 
set of equations for the ft subregion; 


where each of the four columns of [vO)] are the unit 
local deformations associated with the four cases of 
{£>}. Thus, the calculated functions for the first case 
represent the local deformations as a result of ap- 
plied unit axial strain (e) with dimensional units of 
length per unit axial strain. Similarly, the second 
and third cases define the local deformation associ- 
ated with applied bending curvatures (vs, vtj) with 
dimensional units of length per unit bending curva- 
ture. Finally the fourth case describes the local de- 
formation from applied twist rate (6) with dimen- 
sional units of length per unit twist rate. Similarly, 
the stress components of the ft subregion can be 
expressed in terms of a set of unit stresses and [b] 
by substituting Eqns. (15) and (9) into (2. a) 


K (!) 
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where 


where the stiffness matrix is defined as 
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Relations could also be easily developed for the 
displacements (u('), v('), w(')) in terms of a matrix 
of unit displacements and {b}. 

A two-dimensional finite element program was 
written where the cross section is discretized using 
8-node quadrilateral elements and 6-node triangular 
elements. The cross section is defined in a local 
element coordinate system (£ e , rj e ) that can be 
arbitrarily positioned relative to the initial twist ( z ) 
axis of the curvilinear reference frame (g.rj.z) using 
offsets ($/, t ) /). The discretization of thin rectangu- 
lar cross-section (c/t=10) using 40 quadrilateral ele- 
ments is presented in Fig. 2. 


behavior, but not cross section deformation. For 
example, to place a beam in pure torsion with no 
axial strain or bending (0 * 0, e = = 0), one 

must apply an axial force and a general moment with 
bending moment components that satisfy 

P - * 14 ¥l m *34 ^!=*24. (19. a-c) 

M z 4 M z k 4 4 M z /C44 

Positive (negative) ratios of (19. a) are associated 
with the application of an extension force to keep 
the beam from contracting (extending), whereas 
nonzero terms of (19.b,c^ signify that the general 
moment is acting about a vector that is skew (not 
perpendicular) to the cross section. Similarly, to 
place a beam in pure extension with no bending or 
twisting (e *■ 0, 9 = v£ = K-q = 0), one must also apply 
a general moment with components that satisfy 


Behavior of the Pretwisted Beam 


The general behavior of the pretwisted beam 
can be studied by making use of the calculated 
stress distributions (Eq. (16)) and the equilibrium 
equations of the cross section (Eqns. (3.a-d)). 
Thus, 
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(18) 


where the matrix [k] is symmetric based upon re- 
ciprocity. For an untwisted isotropic beam, the 
twist rate is independent of axial strain and bending 
curvatures (^14=^24=^34=0), the diagonal terms 
equal the nominal extensional stiffness (k-\-\ = EA 0 ), 
bending stiffnesses (k22=^no. k33=E/'|o). and 
torsion stiffness (k44=GJo). The remaining off-di- 
agonal terms, which couple extension and bending, 
are the first and second moments of inertia that re-— 
suit when the local axes (5,77) are not coincident 
with the principal axes of the section (ki3=EA 7?0 , 
kl2=-£A* 0 , k23=-ElqT]o)' These last three terms can 
be used to locate the centroid and principal axes of 
the section. For beams exhibiting generally 
anisotropic behavior as a result of material definition 
or from the presence of initial beam twist, the ma- 
trix relation of Eq. (18) will be fully populated. 


_ *13 1 _ *12 Mi ,*14. (20. a-c) 

P *n P *n p *ti 

These ratios agree with the equations developed by 
Lekhnitskii [19] for untwisted generally anisotropic 
beams. 

The behavior of an "unconstrained" anisotropic 
beam with initial twist can be studied by multiplying 
Eq. (18) by the inverse of [k], which results in a 
flexibility relationship. For example, applying an 
axial force (P) produces extension as well as 
bending and twisting that satisfies the following ra- 
tios 


^£.£b., Mu£ a, fi.fii, (21 .a-c) 

e e ^ e H 

where ajj are the components of the flexibility matrix 

[a] (=[k]' 1 ). Similarly, applying a torsion moment 
[M z ) results in extension and bending ratios of the 
form 


£ = il±, Ii*i24, , (22. a-c) 

9 a 44 9 a 44 Q a 44 

Negative (or positive) ratios of (21. c) and (22. a) 
correspond to untwisting (or further twisting) of 
the beam as a result of an applied extension, and 
contraction (or extension) from an applied twist 
moment, respectively. 


Results and Discussion 


The behavior of a "constrained" nonhomo- 
geneous anisotropic beam with initial twist can be 
studied by using Eq. (18) directly, where forces 
and/or moments are applied to restrict global beam 


Three sets of numerical results are presented to 
illustrate the capabilities of the current analytical 
model and how the interaction of material and 
pretwist definitions effect the extension-twist be- 
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havior of thin-wall composite beams. Initially, a vali- 
dation study is performed. Second, the extension- 
twist coupling behavior of a pretwisted 
graphite/epoxy strip, which is geometrically similar 
to a jet turbo-fan blade, is studied for an asymmetric 
[±0] and symmetric ([±0]$) stacking sequence. 
Finally,' a pretwisted beam with a thin-wall D-section 
composed of a graphite/epoxy woven fabric is ana- 
lyzed for different ply orientations and pretwist 
definitions. This beam is geometrically similar to 
blades used on tilt-rotor aircraft (NASA XV-15). 

Verification Studies 

Since published results do not currently exist for 
the behavior of pretwisted thin-wall advanced 
composite beams, the current model was validated 
using published results for pretwisted isotropic 
beams and flat (untwisted) laminated composite 
plates. Because the current approach will reduce 
identically to the elasticity model of ref. [4] for 
pretwisted isotropic beams, it is not necessary in 
this limited space to show that the current model is 
in exact agreement with the isotropic models of 
Refs. [3-9,12] for helical beams having a thin ellipti- 
cal cross-section. Interested readers should refer to 
[4] for a detailed discussion and numerical results 
that include the variation of section and extension- 
twist coupling properties with pretwist definition 
(initial twist level, initial twist axis location). 

The extension-twist coupling behavior of an 
axially-loaded flat laminated plate (a=0, 0) with 

an asymmetric stacking sequence [ro] was analyzed 
and the results were compared to classical lamina- 
tion plate theory predictions [20], The thin plate 
cross section (c/f=l0) is composed of two plies of 
unidirectional graphite/epoxy fibers, where each ply 
is defined using 20 quadrilateral elements (see Fig. 
2). The properties of the fiber system (Table 1) are 
defined relative to an orthogonal reference frame 
(1,2,3) where the 1-axis is coincident with the fiber 
direction. The 1-2 plane of the fiber system is paral- 
lel to the z-q plane of the plate cross section and the 
3 axis is coincident with the r\ axis. A ply angle (0) 
is used to locate the fiber direction (1-axis) relative 
to the (z) axis, where a positive angle is defined as a 
counterclockwise rotation about the tj axis for the 
upper ply and a counterclockwise rotation about the 
-jj axis for the lower ply. If 0=0, then the fibers of 
both the upper and lower plies are parallel to the z- 
axis 

The nondimensionalized ratio of the plate 
untwist (8) for a given extensional strain (e) as a re- 
sult of an applied force (P) can be expressed (from 
Eq. 21. c) as; 

££«£l4_£_ (23. a) 

e ^ 


A similar relationship can be developed from classical 
laminated plate theory by making use of the inverse 
of the laminate stiffness matrix (commonly called 
the " A-B-D“ matrix); 



(23 b) 

where An* and Bis* are the coefficients from the 
inverse of the laminate stiffness matrix that are as- 
sociated with the amount of extension strain and 
twist curvature from an applied in-plane force (N x ), 
respectively. In Fig. 3, the variation of this 
twist/extension ratio as a function of ply angle (0) is 
presented for the current approach and for classical 
laminated plate theory. These results illustrate that 
the current approach is in exact agreement with the 
classical laminated plate theory and that the maxi- 
mum amount of positive (8> 0) or negative (0<O) 
plate twisting for a minimum amount of extension 
strain (e) occurs when (0) is equal to 12° or -12°, 
respectively. Furthermore, orientating the fibers at 
either 0° or 90° will produce plate extension only 
(9=0). 

Pretwisted Laminated Composite.-Stria 

The extension-twist coupling behavior of a 
pretwisted graphite/epoxy strip, which is geometri- 
cally similar to aviation propeller and jet turbo-fan 
blades (0.1<ac<0.2), is studied for an asymmetric 
[±<z>] and symmetric ([±0]$) stacking sequence. The 
thin 2-ply asymmetric cross-section of the verifica- 
tion study (see Fig. 2) will again be analyzed except 
now the strip has initial twist defined about the 
centroidal axis (£prj^0). In Fig. 4, the nondimen- 
sionalized extension/twist flexibility coefficient 
(a-uEiic 3 ) is presented as a function of ply angle 
(0) and nondimensionalized pretwist ( ac ). For low 
pretwist levels (ac<0.0l), the (ai4Eiic^) curve 
undergoes a slight downward shift. This shift is a 
result of untwisting of the axially loaded pretwisted 
beam and thus the flexibility coefficient becomes 
more negative for negative ply angles, less positive 
for positive ply angles, and the zero coefficient 
value shifts to a positive ply angle. For moderate 
pretwist levels (ac=0.l0), the downward shift of 
the curve is substantial, where the initial twist re- 
lated coupling values are of the same order as the 
material related coupling effects and the ply angles 
for maximum coupling have increased by 4°. 
Furthermore, the ply angles for zero coupling has 
shifted from 0° and -90° to 11° and 60°. Finally, 


for large pretwist levels (ac= 0.20, 0.30), the 
downward shift of the curve can be large enough 
that eventually no positive or zero coupling exists 
and thus an axially loaded pretwisted strip will al- 
ways untwist independent of ply angle definition. 

The variation in the extension stiffness as a 
function of fiber angle and pretwist level is pre- 
sented in Fig. 5, where EAq (=1/ai i 0 ) is defined as 
the extension stiffness of a unidirectional flat strip 
(cr=0= 0). The downward shift of the curves for 
moderate to highly pretwisted strips was expected, 
since it is well known that adding initial twist to an 
isotropic beam will lower the extension stiffness 
(see [4]). The shift of the relative maximum from 
(0=0) to a small positive ply angle is a result of using 
the material coupling to counteract the 
twist/extension coupling associated with beam initial 
twist. The ply angles for maximum extension stiff- 
ness are equal to, for small initial twist levels, the 
zero values of twist/extension flexibility (from Fig. 
4). In Fig. 6, the twist/extension ratio (0c/e) is pre- 
sented, where the maximum untwisting (negative) 
for minimum extension strain generally occurs for 
ply angles between 0° and -12°. The ratios for 
maximum (positive) twisting undergo a significant 
reduction and the ply angle definitions for zero 
coupling are highly dependent upon the initial twist 
level. Finally, for small to moderate pretwist levels, 
the slopes of the curves can be extremely steep and 
thus small ply angle changes will cause large changes 
in the resulting ratios. 

The torsion stiffness is presented in Fig. 7, 
where GJo (=1/a44 0 ) is the nominal torsion stiff- 
ness of a unidirectional flat strip (a=0=O). For strips 
with small initial twist (ac<0.01), the variation in the 
torsion stiffness is nearly identical to a flat strip, 
where the relative maximums occur at <$>=±45°. For 
moderate to large pretwist levels, the torsion stiff- 
ness will increase (as expected), where the largest 
increase occurs for small ply angles (near zero) with 
the magnitude of the increase tapering off near 
0=±45°. These results are most interesting in that 
ply angles less than ±45° undergo the largest per- 
centage increase. 

A contour plot of the twist/extension ratio 
(©c/e) as a function of ply angle (<p) and initial twist 
axis offset along the chord (£/*0, rj/= 0) for the 2-ply 
asymmetric strip with initial twist (ac=0.l0) is pre- 
sented in Fig. 8. The dashed contour lines represent 
a mid-level between the solid contour lines. The 
maximum negative ratios (-35) occur when qi=0 and 

0=-5°, whereas a local maximum positive ratio (+15) 
can be found at %/=0 and <£=20° but a global positive 
maximum will occur at large values of |/c. Thus, an 
initially negative ratio can always be made zero or 
positive by shifting the initial twist axis outward 


from the centroid. Finally, a region exists (<£=45°) 
where the twist/extension ratio is independent of 
initial twist axis location. 

This thin rectangular graphite/epoxy cross- 
section (c/t=10) was also studied using a 4-ply sym- 
metric ([±0)s) stacking sequence, where each ply is 
defined using 20 quadrilateral elements (80 elements 
total). and a positive ply angle {<t>) on the top and 
bottom plies is defined as a counterclockwise rota- 
tion about the rj axis and on the two inner plies as a 
counterclockwise rotation about the -77 axis. Initially, 
the initial twist is defined to act through the section 
centroid (|/=tjt= 0). The twist/extension ratio (9c/e) 
is presented in Fig. 9, where extension-torsion be- 
havior occurs as a result of initial twist, since a flat 
strip with a symmetric ply lay-up will always have 
zero coupling. It is apparent that maximum untwist- 
ing occurs with a 0° ply angle and for nonzero ply 
configurations the amount of untwist is greatly re- 
duced. Comparing Figs. 6 and 9, both cross-sec- 
tions have the same magnitudes for 0=0°, but the 
asymmetric ply definition provides greater freedom 
for controlling (optimizing) twist/extension cou- 
pling. 

The variation in the extension (EA/EA 0 ) and 
torsion stiffnesses (GJ/GJ 0 ) are presented in Figs. 
10 and 11, respectively, where the curves are sym- 
metric with respect to 0=0°. Comparing the varia- 
tion of the extension stiffness with the 2-ply asym- 
metric results (Fig. 5), it is readily apparent that 
both results agree for 0= 0° and the symmetric lay- 
up results are nearly constant for small ply angles (- 
10°<<?<100), whereas the asymmetric results are 
extremely sensitive to small ply angles. Similarly, the 
increase in the torsion stiffness agrees with the 
asymmetric ply lay-up results (Fig. 6) for 0=0°, 
0>45°, and o<-45°, but the asymmetric ply lay-up 
can have a much greater torsion stiffness for 
O°<0<45°. 

A contour plot of the twist/extension ratio 
(6c/e) as a function of ply angle (0) and initial twist 
axis offset along the chord (£,'* 0 , 77 /=0) with initial 
twist (ac=0.i0) is presented in Fig. 12. Again, the 
dashed contour lines represent a mid-level between 
the solid contour lines. The maximum negative (-32) 
extension/twist ratios occur with £/= 0 and 0=0°, 
where changing the ply angle or initial twist axis lo- 
cation will only increase the ratio. The values for 
zero coupling are independent of ply angle and thus 
negative coupling will always exist as long as (- 
0.28<<-/c <0.28), but this coupling can be small for 
large ply angles. 

Pretwisted Laminated Composite D-Section 


The extension-twist coupling behavior of a 
pretwisted beam with a thin-wall D-section com- 
posed of a constant thickness graphite/epoxy wo- 
ven fabric is analyzed for different ply [0] orienta- 
tions and pretwist definitions. This beam is geo- 
metrically similar to the structural box beam used in 
tilt-rotor aircraft blades (0.01<ac<0.i0). The D- 
section geometry is presented in Fig. 13, where 150 
quadrilateral elements are used to define the single- 
cell planform. The material properties of the 
graphite/epoxy woven cloth are given in Table 1 and 
a positive ply angle <p is defined as a counterclock- 
wise rotation about a outward normal (n) that is 
perpendicular section surface. Thus, the ply angle is 
defined relative to the local element coordinate sys- 
tem to simulate the wrapping of the woven cloth. 
Even though the wall laminate properties are sym- 
metric, as a result of treating the woven cloth as a 
single ply, the resulting beam will experience 
twist/extension coupling because the effective 
cross-section has an asymmetric material definition. 

The twist/extension ratio (6c/e) is presented in 
Fig. 14, where the complete behavioral range can be 
described over a 90° ply angle period, instead of 
the 180° period typical of uniaxial fibers. The flat 
(a=0) beam section obviously behaves in a fashion 
similar to the asymmetric ply definition (Fig. 6) with 
positive and negative coupling for positive (O<0 
<45°) and negative (-45 °<<p <0) ply angles, respec- 
tively. Adding initial twist, about the cross-section 
centroid (Sf=0.523c, tj/= 0), shifts the curves down- 
ward with the largest changes in the magnitude oc- 
curring with (-15°<0 <15°). Finally, the magnitudes 
are significantly less than the magnitudes associated 
with the uniaxial fiberous material (asymmetric; Fig. 
6, symmetric; Fig. 9) because the woven material 
has a lower extension stiffness to shear stiffness 
(E/G) ratio. A contour plot of the twist/extension 
ratio (0c/e) as a function of ply angle (0) and initial 
twist axis offset (0<|/c<l, 77 7=0) for (ac=0.l0) is 
presented in Fig. 15. This plot has many of the 
same features of Fig. 8, with maximum positive and 
negative regions occurring with the initial twist axis 
acting through the cross-section centroid. The zero 
coupling, which encircles the negative coupling, can 
be approximated as a series of parallel (vertical) 
lines, as opposed to the two horizontal parallel lines 
for the symmetric strip (Fig. 12). Thus, the zero 
coupling is nearly independent of the initial twist 
axis location for most the chord length. 
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Table 1: Material properties of T300/5208 

graphite/epoxy unidirectional fibers and 
woven cloth. 

unidirectional fibers woven cloth 


£11 

132.2 GPa 

80.32 GPa 

E 22 

10.75 GPa 

80.32 GPa 

^33 

10.75 GPa 

10.75 GPa 

Gl 2 ■ ^13 = G 23 

5.65 GPa 

5.65 GPa 

vi 2 

0.2M 

0.050 

vi 3 

0.239 

0.239 

v 23 

0.400 

0.239 



Fig. 1 Nonhomogeneous anisotropic beam with 
initial twist. 



Fig. 2 Finite element discretization of a thin 
rectangular cross section (c/t=i0). 
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Fig. 3 Twist/extension ratio for a 2-ply 

asymmetric graphite/epoxy untwisted 
strip. 
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Fig. 4 Nondimensionalized extension/twist 
coefficient of a 2-ply asymmetric [±0] 
graphite/epoxy strip with initial twist about 
the centroidal axis (£/= rjf=0). 


Fig. 6 Twist/extension ratio for a 2-ply 

asymmetric [±0] graphite/epoxy strip with 
initial twist about the centroidal axis 

(&=W=0). 



Fig. 5 Axial stiffness of a 2-ply asymmetric [±0] Fig. 7 Torsion stiffness of a 2-ply asymmetric [±0] 
graphite/epoxy strip with initial twist about graphite/epoxy strip with initial twist about 

the centroidal axis (€/=rj/=0). the centroidal axis (^rj^O). 
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Fig. 8 Contour plot of the twist/extension ratio 
for a 2-ply asymmetric [±<j>] graphite/epoxy 
pretwisted strip (ac=0.10). 



Fig. 9 Twist/extension ratio for a 4-ply symmetric 
[±0] s graphite/epoxy strip with initial twist 
about the centroidal axis (£/=tj/= 0). 



Fig. 10 Axial stiffness of a 4-ply symmetric [±$>] s 

graphite/epoxy strip with initial twist about 
the centroidal axis (?/=rjp0). 



Ply Angle ♦ {degree) 

Fig. 1 1 Torsion stiffness of a 4-ply symmetric [±$>]s 
graphite/epoxy strip with initial twist about 
the centroidal axis (|/=rj/=0). 
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Fig. 12 Contour plot of twist/extension ratio for a 
4-ply symmetric [±0] s graphite/epoxy 
pretwisted strip (ac=0.l0). 


Fig. 14 Twist/extension ratio for a woven 

graphite/epoxy D-section with initial twist 
about the centroidal axis (£/=0.523c, rj/= 0). 


, initial twist axis 


Fig. 13 Thin wall D-section composed of a single 
layer [ 0 ] of graphite/epoxy woven cloth. 
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Fig. 15 Contour plot of twist/extension ratio for a 
woven graphite/epoxy D-section with 
(ac=0.10). 








